From 4df1bc87c94f9b740d814c417462443cbe6ebf85 Mon Sep 17 00:00:00 2001 From: Miguel de la Varga Date: Mon, 27 Jul 2026 14:19:17 +0200 Subject: [PATCH] [ENH] Introduce `FiniteFault` model with validation, serialization, and spline taper support Adds a Pydantic-based `FiniteFault` dataclass to define finite faults with validated, JSON-serializable fields. Supports cubic, quadratic, and spline tapers, with validation for radii, rotation, and spline control points. Refactors slip calculation logic to use `FiniteFault` and centralizes controls for consistency. Includes tests for serialization, validation, and taper configurations, along with updated documentation. --- docs/implementation/finite_faults.md | 135 ++++++++++++++++++ gempy_engine/core/data/__init__.py | 1 + gempy_engine/core/data/finite_fault.py | 86 +++++++++++ gempy_engine/modules/faults/finite_faults.py | 97 ++++--------- .../test_faults/test_finite_fault_api.py | 2 +- .../test_finite_fault_integration.py | 6 +- .../test_core/test_finite_fault.py | 81 +++++++++++ 7 files changed, 338 insertions(+), 70 deletions(-) create mode 100644 docs/implementation/finite_faults.md create mode 100644 gempy_engine/core/data/finite_fault.py create mode 100644 tests/test_common/test_core/test_finite_fault.py diff --git a/docs/implementation/finite_faults.md b/docs/implementation/finite_faults.md new file mode 100644 index 00000000..8f2b5cb4 --- /dev/null +++ b/docs/implementation/finite_faults.md @@ -0,0 +1,135 @@ +# Finite Faults Implementation + +## Goal + +Define finite faults as validated, serializable input data and use that +definition to taper fault displacement on the interpolated fault surface. + +The implementation is based on the projection and UV taper prototype in +`gempy_engine/modules/faults/finite_faults.py`. The prototype is useful, but it +is not yet connected to the engine's older callable-based finite-fault hook. + +## Status + +### Phase 1: Serializable input definition + +- [x] Introduce a Pydantic dataclass for the finite-fault definition. +- [x] Keep persisted fields JSON-native rather than storing NumPy arrays. +- [x] Serialize taper types as stable string values. +- [x] Validate centers, radii, rotation, and spline control points. +- [x] Add JSON and Python round-trip tests using `pydantic.TypeAdapter`. +- [x] Preserve the prototype import path and numerical API. + +### Phase 2: Projection correctness + +- [ ] Remove the incorrect fixed half-step from plane projection. +- [ ] Define explicit behavior for near-zero gradients. +- [ ] Test exact projection onto a plane and projection idempotency. +- [ ] Decide whether nonlinear fields need iterative re-evaluation. +- [ ] Fix the dense-grid gradient accessor before using it for projection. + +### Phase 3: Stack wiring + +- [ ] Replace the callable in `FiniteFaultData` with the declarative definition. +- [ ] Associate at most one finite-fault definition with each fault stack. +- [ ] Validate that definitions are attached only to `StackRelationType.FAULT` stacks. +- [ ] Make scalar gradients available when a finite-fault stack is evaluated. +- [ ] Pass the fault scalar field, gradients, and surface isovalue to the taper operation. +- [ ] Apply the taper to the fault drift before dependent stacks are interpolated. +- [ ] Cover both sequential and flat-stack interpolation paths. + +### Phase 4: Integration and backend support + +- [ ] Add a numerical integration test for a dependent stratigraphic stack. +- [ ] Verify that displacement reaches zero at the finite-fault tips. +- [ ] Document the approximately planar local-frame limitation. +- [ ] Define and test NumPy and PyTorch backend behavior. +- [ ] Add the finite-fault definition to the server payload when stack data is exposed there. + +## Input Contract + +`FiniteFault` is a frozen Pydantic dataclass. Its persisted representation +contains only JSON-compatible values: + +| Field | Type | Meaning | +| --- | --- | --- | +| `center` | 3-tuple of floats | Point at the center of the finite-fault footprint | +| `strike_radius` | positive float or 2-tuple | Positive and negative strike radii | +| `dip_radius` | positive float or 2-tuple | Positive and negative dip radii | +| `taper` | `cubic`, `quadratic`, or `spline` | Slip taper profile | +| `rotation_deg` | float | In-plane rotation in degrees | +| `spline_control_points` | optional sequence of 2-tuples | Distance-to-slip profile for a spline taper | + +For an asymmetric radius, tuple order is `(positive_direction, +negative_direction)`. All radii must be finite and greater than zero. + +Spline control points use `(normalized_distance, slip_multiplier)`. Distances +must be strictly increasing from `0` to `1`, and multipliers must remain in the +range `[0, 1]`. Omitting the points for a spline taper selects the engine's +default profile. Supplying spline points for another taper is invalid. + +`normal_radius` is intentionally not part of this contract. The UV workflow +projects points onto the fault surface, so normal distance is not part of its +two-dimensional footprint. A volumetric ellipsoid would be a separate model. + +## Serialization + +Pydantic dataclasses use `TypeAdapter` for serialization and deserialization: + +```python +from pydantic import TypeAdapter + +from gempy_engine.core.data.finite_fault import FiniteFault, TaperType + +adapter = TypeAdapter(FiniteFault) + +finite_fault = FiniteFault( + center=(0.0, 0.0, 0.0), + strike_radius=(2.0, 1.0), + dip_radius=0.75, + taper=TaperType.SPLINE, + rotation_deg=15.0, + spline_control_points=( + (0.0, 1.0), + (0.5, 0.8), + (1.0, 0.0), + ), +) + +payload: bytes = adapter.dump_json(finite_fault) +restored: FiniteFault = adapter.validate_json(payload) +``` + +Equivalent JSON: + +```json +{ + "center": [0.0, 0.0, 0.0], + "strike_radius": [2.0, 1.0], + "dip_radius": 0.75, + "taper": "spline", + "rotation_deg": 15.0, + "spline_control_points": [ + [0.0, 1.0], + [0.5, 0.8], + [1.0, 0.0] + ] +} +``` + +## Known Prototype Issues + +- `project_points_onto_surface` currently moves points only halfway to a linear plane. +- Near-zero gradients are silently replaced with a denominator of one. +- `FiniteFault.calculate_slip` expects callers to project points separately. +- The local strike/dip frame is constant and therefore approximates curved faults. +- The engine-integrated `FiniteFaultData` loses its callable when serialized. +- `ScalarFieldOutput.exported_fields_dense_grid` currently returns scalar values as gradients. +- Existing integration assertions do not verify the projected surface residual or expected geometry. + +## Design Decisions + +- The finite-fault definition is declarative; serialized callables are not supported. +- NumPy conversion happens at the numerical boundary, not in persisted fields. +- Geometry belongs to the fault stack that defines it, not each destination stack affected by it. +- Projection algorithm settings are evaluation concerns and are not part of geological input data. diff --git a/gempy_engine/core/data/__init__.py b/gempy_engine/core/data/__init__.py index 1b90dd67..3bfa6adc 100644 --- a/gempy_engine/core/data/__init__.py +++ b/gempy_engine/core/data/__init__.py @@ -1,3 +1,4 @@ +from .finite_fault import FiniteFault, TaperType from .tensors_structure import TensorsStructure from .kernel_classes.orientations import Orientations, OrientationsInternals from .kernel_classes.surface_points import SurfacePoints, SurfacePointsInternals diff --git a/gempy_engine/core/data/finite_fault.py b/gempy_engine/core/data/finite_fault.py new file mode 100644 index 00000000..e67a884c --- /dev/null +++ b/gempy_engine/core/data/finite_fault.py @@ -0,0 +1,86 @@ +from __future__ import annotations + +import math +from enum import Enum +from typing import Any + +import numpy as np +from pydantic import ConfigDict, field_validator, model_validator +from pydantic.dataclasses import dataclass + + +class TaperType(str, Enum): + CUBIC = "cubic" + QUADRATIC = "quadratic" + SPLINE = "spline" + + +Radius = float | tuple[float, float] +SplineControlPoints = tuple[tuple[float, float], ...] + + +@dataclass(frozen=True, config=ConfigDict(extra="forbid")) +class FiniteFault: + """Serializable definition of an approximately planar finite fault.""" + + center: tuple[float, float, float] + strike_radius: Radius = 1.0 + dip_radius: Radius = 1.0 + taper: TaperType = TaperType.CUBIC + rotation_deg: float = 0.0 + spline_control_points: SplineControlPoints | None = None + + @field_validator("center", "strike_radius", "dip_radius", "spline_control_points", mode="before") + @classmethod + def _convert_numpy_arrays(cls, value: Any) -> Any: + return value.tolist() if isinstance(value, np.ndarray) else value + + @field_validator("center") + @classmethod + def _validate_center(cls, value: tuple[float, float, float]) -> tuple[float, float, float]: + if not all(math.isfinite(component) for component in value): + raise ValueError("center coordinates must be finite") + return value + + @field_validator("strike_radius", "dip_radius") + @classmethod + def _validate_radius(cls, value: Radius) -> Radius: + radii = value if isinstance(value, tuple) else (value,) + if not all(math.isfinite(radius) and radius > 0 for radius in radii): + raise ValueError("radii must be finite and greater than zero") + return value + + @field_validator("rotation_deg") + @classmethod + def _validate_rotation(cls, value: float) -> float: + if not math.isfinite(value): + raise ValueError("rotation_deg must be finite") + return value + + @model_validator(mode="after") + def _validate_spline_control_points(self) -> "FiniteFault": + points = self.spline_control_points + if points is None: + return self + if self.taper is not TaperType.SPLINE: + raise ValueError("spline_control_points require a spline taper") + if len(points) < 2: + raise ValueError("spline_control_points require at least two points") + + distances = tuple(point[0] for point in points) + multipliers = tuple(point[1] for point in points) + if distances[0] != 0.0 or distances[-1] != 1.0: + raise ValueError("spline distances must start at 0 and end at 1") + if not all(math.isfinite(value) for point in points for value in point): + raise ValueError("spline control points must be finite") + if not all(left < right for left, right in zip(distances, distances[1:])): + raise ValueError("spline distances must be strictly increasing") + if not all(0.0 <= multiplier <= 1.0 for multiplier in multipliers): + raise ValueError("spline multipliers must be between 0 and 1") + return self + + def calculate_slip(self, points: np.ndarray, normal: np.ndarray) -> np.ndarray: + """Calculate the prototype slip multiplier for the supplied points.""" + from gempy_engine.modules.faults.finite_faults import calculate_slip + + return calculate_slip(self, points=points, normal=normal) diff --git a/gempy_engine/modules/faults/finite_faults.py b/gempy_engine/modules/faults/finite_faults.py index ba65cc53..717b6950 100644 --- a/gempy_engine/modules/faults/finite_faults.py +++ b/gempy_engine/modules/faults/finite_faults.py @@ -1,73 +1,38 @@ import numpy as np -from enum import Enum -from dataclasses import dataclass -from typing import Optional, Union, Tuple from scipy.interpolate import make_interp_spline - -class TaperType(Enum): - CUBIC = "cubic" - QUADRATIC = "quadratic" - SPLINE = "spline" - - -@dataclass -class FiniteFault: - """ - Elegant API for defining finite faults. - - Args: - center: Center of the fault in 3D space. - strike_radius: Radius along the strike direction (u). - Can be a single float or a tuple (positive_u, negative_u) for anisotropy. - dip_radius: Radius along the dip direction (v). - Can be a single float or a tuple (positive_v, negative_v) for anisotropy. - normal_radius: Radius along the normal direction (w). Default 1.0. - taper: The tapering function to use. - spline_control_points: If taper is SPLINE, these points define the curve. - """ - center: np.ndarray - strike_radius: Union[float, Tuple[float, float]] = 1.0 - dip_radius: Union[float, Tuple[float, float]] = 1.0 - normal_radius: Optional[Union[float, Tuple[float, float]]] = None - taper: TaperType = TaperType.CUBIC - rotation: float = 0.0 - spline_control_points: Optional[np.ndarray] = None - - def __post_init__(self): - if self.taper == TaperType.SPLINE and self.spline_control_points is None: - # Default bell-shaped spline if none provided - self.spline_control_points = np.array([ - [0.0, 1.0], - [0.2, 0.95], - [0.5, 0.5], - [0.8, 0.05], - [1.0, 0.0] - ]) - - def calculate_slip(self, points: np.ndarray, normal: np.ndarray) -> np.ndarray: - """ - High-level method to calculate slip multiplier for given points. - """ - u, v, w = get_local_frame(normal, angle_deg=self.rotation) - d = get_ellipsoid_distance( - points=points, - center=self.center, - u=u, v=v, w=w if self.normal_radius is not None else None, - a=self.strike_radius, - b=self.dip_radius, - c=self.normal_radius if self.normal_radius is not None else 1.0 - ) - - if self.taper == TaperType.CUBIC: - return cubic_hermite_taper(d) - elif self.taper == TaperType.QUADRATIC: - return quadratic_taper(d) - elif self.taper == TaperType.SPLINE: - return spline_taper(d, self.spline_control_points) - else: - raise ValueError(f"Unknown taper type: {self.taper}") +from gempy_engine.core.data.finite_fault import FiniteFault, TaperType + +_DEFAULT_SPLINE_CONTROL_POINTS = np.array([ + [0.0, 1.0], + [0.2, 0.95], + [0.5, 0.5], + [0.8, 0.05], + [1.0, 0.0] +]) + + +def calculate_slip(finite_fault: FiniteFault, points: np.ndarray, normal: np.ndarray) -> np.ndarray: + """Calculate a slip multiplier using the prototype local-plane model.""" + u, v, _ = get_local_frame(normal, angle_deg=finite_fault.rotation_deg) + distance = get_ellipsoid_distance( + points=points, + center=np.asarray(finite_fault.center), + u=u, + v=v, + a=finite_fault.strike_radius, + b=finite_fault.dip_radius, + ) + + if finite_fault.taper is TaperType.CUBIC: + return cubic_hermite_taper(distance) + if finite_fault.taper is TaperType.QUADRATIC: + return quadratic_taper(distance) + if finite_fault.taper is TaperType.SPLINE: + control_points = finite_fault.spline_control_points or _DEFAULT_SPLINE_CONTROL_POINTS + return spline_taper(distance, np.asarray(control_points)) + raise ValueError(f"Unknown taper type: {finite_fault.taper}") def get_local_frame(normal: np.ndarray, angle_deg: float = 0) -> tuple[np.ndarray, np.ndarray, np.ndarray]: diff --git a/tests/test_common/test_api/test_faults/test_finite_fault_api.py b/tests/test_common/test_api/test_faults/test_finite_fault_api.py index c969c986..9219dd53 100644 --- a/tests/test_common/test_api/test_faults/test_finite_fault_api.py +++ b/tests/test_common/test_api/test_faults/test_finite_fault_api.py @@ -84,7 +84,7 @@ def test_finite_fault_rotation(): center=center, strike_radius=2.0, dip_radius=1.0, - rotation=90 # Rotate 90 deg + rotation_deg=90 # Rotate 90 deg ) # Normal along Z => default strike u=[1,0,0], dip v=[0,1,0] diff --git a/tests/test_common/test_api/test_faults/test_finite_fault_integration.py b/tests/test_common/test_api/test_faults/test_finite_fault_integration.py index 986f1019..9b94e10a 100644 --- a/tests/test_common/test_api/test_faults/test_finite_fault_integration.py +++ b/tests/test_common/test_api/test_faults/test_finite_fault_integration.py @@ -68,9 +68,9 @@ def test_finite_fault_integration_suite(finite_fault_setup): FiniteFault(center=center, strike_radius=0.8, dip_radius=0.4, taper=TaperType.CUBIC), FiniteFault(center=center, strike_radius=(1.0, 0.5), dip_radius=0.6, taper=TaperType.QUADRATIC), FiniteFault(center=center, strike_radius=0.7, dip_radius=0.7, taper=TaperType.SPLINE, spline_control_points=cp_bell), - FiniteFault(center=center, strike_radius=0.8, dip_radius=0.4, taper=TaperType.CUBIC, rotation=45), - FiniteFault(center=center, strike_radius=0.8, dip_radius=1, taper=TaperType.SPLINE, spline_control_points=cp_plateau, rotation=30), - FiniteFault(center=center, strike_radius=(1.2,0.8), dip_radius=(2,1), taper=TaperType.SPLINE, spline_control_points=cp_plateau, rotation=30) + FiniteFault(center=center, strike_radius=0.8, dip_radius=0.4, taper=TaperType.CUBIC, rotation_deg=45), + FiniteFault(center=center, strike_radius=0.8, dip_radius=1, taper=TaperType.SPLINE, spline_control_points=cp_plateau, rotation_deg=30), + FiniteFault(center=center, strike_radius=(1.2,0.8), dip_radius=(2,1), taper=TaperType.SPLINE, spline_control_points=cp_plateau, rotation_deg=30) ] if not plot_pyvista: diff --git a/tests/test_common/test_core/test_finite_fault.py b/tests/test_common/test_core/test_finite_fault.py new file mode 100644 index 00000000..f909d908 --- /dev/null +++ b/tests/test_common/test_core/test_finite_fault.py @@ -0,0 +1,81 @@ +import json + +import numpy as np +import pytest +from pydantic import TypeAdapter, ValidationError + +from gempy_engine.core.data.finite_fault import FiniteFault, TaperType + + +def test_finite_fault_json_round_trip(): + finite_fault = FiniteFault( + center=np.array([1.0, 2.0, 3.0]), + strike_radius=(2.0, 1.0), + dip_radius=0.75, + taper=TaperType.SPLINE, + rotation_deg=15.0, + spline_control_points=np.array([ + [0.0, 1.0], + [0.5, 0.8], + [1.0, 0.0], + ]), + ) + adapter = TypeAdapter(FiniteFault) + + payload = adapter.dump_json(finite_fault) + restored = adapter.validate_json(payload) + + assert restored == finite_fault + assert json.loads(payload) == { + "center" : [1.0, 2.0, 3.0], + "strike_radius" : [2.0, 1.0], + "dip_radius" : 0.75, + "taper" : "spline", + "rotation_deg" : 15.0, + "spline_control_points": [[0.0, 1.0], [0.5, 0.8], [1.0, 0.0]], + } + + +@pytest.mark.parametrize( + "field,value", + [ + ("strike_radius", 0.0), + ("strike_radius", (1.0, -1.0)), + ("dip_radius", float("inf")), + ], +) +def test_finite_fault_rejects_invalid_radii(field, value): + with pytest.raises(ValidationError, match="radii must be finite and greater than zero"): + FiniteFault(center=(0.0, 0.0, 0.0), **{field: value}) + + +@pytest.mark.parametrize( + "points,error", + [ + (((0.1, 1.0), (1.0, 0.0)), "must start at 0"), + (((0.0, 1.0), (0.5, 0.8), (0.5, 0.4), (1.0, 0.0)), "strictly increasing"), + (((0.0, 1.1), (1.0, 0.0)), "must be between 0 and 1"), + ], +) +def test_finite_fault_rejects_invalid_spline_control_points(points, error): + with pytest.raises(ValidationError, match=error): + FiniteFault( + center=(0.0, 0.0, 0.0), + taper=TaperType.SPLINE, + spline_control_points=points, + ) + + +def test_finite_fault_rejects_spline_points_for_polynomial_taper(): + with pytest.raises(ValidationError, match="require a spline taper"): + FiniteFault( + center=(0.0, 0.0, 0.0), + spline_control_points=((0.0, 1.0), (1.0, 0.0)), + ) + + +def test_finite_fault_rejects_unknown_fields_on_deserialization(): + adapter = TypeAdapter(FiniteFault) + + with pytest.raises(ValidationError, match="Unexpected keyword argument"): + adapter.validate_python({"center": [0.0, 0.0, 0.0], "unknown": True})