From 8c26203a400ab26ff613e312b7774300ade64a9a Mon Sep 17 00:00:00 2001 From: Miguel de la Varga Date: Mon, 27 Jul 2026 14:52:02 +0200 Subject: [PATCH] [ENH] Add finite-fault support to server schema, validation, and PyTorch backend Extends server schema to support finite-fault definitions with validation for fault relations and input data, including JSON serialization. Integrates finite-fault handling into PyTorch backend with device consistency and added tests. Updates documentation with server payload structure and backend behavior details. --- docs/implementation/finite_faults.md | 75 +++++++++++++++++-- .../core/data/input_data_descriptor.py | 15 +++- .../kernel_classes/server/input_parser.py | 9 +++ .../test_finite_fault_stack_wiring.py | 50 ++++++++++++- tests/test_server/test_finite_fault_input.py | 50 +++++++++++++ 5 files changed, 190 insertions(+), 9 deletions(-) create mode 100644 tests/test_server/test_finite_fault_input.py diff --git a/docs/implementation/finite_faults.md b/docs/implementation/finite_faults.md index eed809f5..0fd8c9e0 100644 --- a/docs/implementation/finite_faults.md +++ b/docs/implementation/finite_faults.md @@ -41,11 +41,11 @@ for Phase 4. ### 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. +- [x] Add a numerical integration test for a dependent stratigraphic stack. +- [x] Verify that displacement reaches zero at the finite-fault tips. +- [x] Document the approximately planar local-frame limitation. +- [x] Define and test NumPy and PyTorch backend behavior. +- [x] Add the finite-fault definition to the server payload. ## Input Contract @@ -172,11 +172,74 @@ use the non-stacked symbolic evaluator. Ordinary stacks continue to use the optimized stacked evaluator. NumPy and PyKeOps reductions can differ slightly because their floating-point reduction orders are different. +## Server Payload + +The server descriptor accepts the same finite-fault definition under an entry +aligned with the fault stack. The fault dependency matrix is required when +other stacks consume that fault: + +```json +{ + "input_data_descriptor": { + "number_of_points_per_surface": [9, 12], + "number_of_points_per_stack": [9, 12], + "number_of_orientations_per_stack": [1, 2], + "number_of_surfaces_per_stack": [1, 1], + "masking_descriptor": [3, 1], + "faults_relations": [ + [false, true], + [false, false] + ], + "faults_input_data": [ + { + "thickness": null, + "finite_fault": { + "center": [0.0, 0.0, 0.0], + "strike_radius": [2.0, 1.0], + "dip_radius": 0.75, + "taper": "quadratic", + "rotation_deg": 15.0, + "spline_control_points": null + } + }, + null + ] + } +} +``` + +`masking_descriptor` uses the existing numeric `StackRelationType` values; +`3` denotes `FAULT` and `1` denotes `ERODE`. + +## Geometry Limitation + +The current implementation creates one strike/dip frame for the whole finite +fault. Its normal is taken from the valid projected evaluation point nearest the +configured center. This is appropriate for planar and approximately planar +faults. On a strongly curved fault, one fixed frame distorts distances and does +not constitute a true surface UV parameterization. Supporting those faults will +require a spatially varying frame or a mesh-based UV map. + +## Backend Behavior + +Projection and taper evaluation use NumPy and SciPy. NumPy inputs remain NumPy. +PyTorch CPU or GPU tensors are detached and transferred to CPU for geometric +evaluation, and the resulting multiplier is converted back to the active +backend, device, and dtype before it is applied to the fault drift. Flat +NumPy/PyKeOps evaluation is also supported through an isolated, non-stacked +finite-fault evaluation path. + +The CPU conversion means projection and taper geometry are not differentiable +with respect to the fault scalar field, gradients, center, radii, rotation, or +spline points. A differentiable finite-fault workflow would require backend- +native projection and taper functions. NumPy and PyTorch CPU behavior are covered +by tests; GPU placement follows the same conversion path and is exercised only +when CUDA is available. + ## Known Prototype Issues - Standalone `FiniteFault.calculate_slip` expects callers to project points separately; engine stack wiring performs that projection. - The local strike/dip frame is constant and therefore approximates curved faults. -- Existing integration assertions do not verify the projected surface residual or expected geometry. ## Design Decisions diff --git a/gempy_engine/core/data/input_data_descriptor.py b/gempy_engine/core/data/input_data_descriptor.py index 85dbd8a1..d95c2d44 100644 --- a/gempy_engine/core/data/input_data_descriptor.py +++ b/gempy_engine/core/data/input_data_descriptor.py @@ -51,11 +51,24 @@ def from_schema(cls, schema: InputDataDescriptorSchema): # Convert list of ints into list of StackRelationType list_relations: list[StackRelationType] = [StackRelationType(x) for x in schema.masking_descriptor] + faults_relations = None if schema.faults_relations is None else np.array(schema.faults_relations, dtype=bool) + faults_input_data = None + if schema.faults_input_data is not None: + faults_input_data = [ + None if fault_data is None else FaultsData.from_user_input( + thickness=fault_data.thickness, + finite_fault=fault_data.finite_fault, + ) + for fault_data in schema.faults_input_data + ] + stack_structure = StacksStructure( number_of_points_per_stack=np.array(schema.number_of_points_per_stack), number_of_orientations_per_stack=np.array(schema.number_of_orientations_per_stack), number_of_surfaces_per_stack=np.array(schema.number_of_surfaces_per_stack), - masking_descriptor=list_relations + masking_descriptor=list_relations, + faults_relations=faults_relations, + faults_input_data=faults_input_data, ) return cls(tensors_structure=tensor_structure, stack_structure=stack_structure) diff --git a/gempy_engine/core/data/kernel_classes/server/input_parser.py b/gempy_engine/core/data/kernel_classes/server/input_parser.py index 221faf74..659af4bf 100644 --- a/gempy_engine/core/data/kernel_classes/server/input_parser.py +++ b/gempy_engine/core/data/kernel_classes/server/input_parser.py @@ -1,6 +1,8 @@ # Check if pydantic is installed and import it from typing import Optional +from ...finite_fault import FiniteFault + try: from pydantic import BaseModel, Field except ImportError: @@ -27,12 +29,19 @@ class InterpolationInputSchema(BaseModel): grid: GridSchema +class FaultsDataSchema(BaseModel): + thickness: Optional[float] = None + finite_fault: Optional[FiniteFault] = None + + class InputDataDescriptorSchema(BaseModel): number_of_points_per_surface: list[int] number_of_points_per_stack: list[int] number_of_orientations_per_stack: list[int] number_of_surfaces_per_stack: list[int] masking_descriptor: list[int] # * StackRelationType + faults_relations: Optional[list[list[bool]]] = None + faults_input_data: Optional[list[Optional[FaultsDataSchema]]] = None class GemPyInput(BaseModel): diff --git a/tests/test_common/test_api/test_faults/test_finite_fault_stack_wiring.py b/tests/test_common/test_api/test_faults/test_finite_fault_stack_wiring.py index 785e1cb1..8c0bab20 100644 --- a/tests/test_common/test_api/test_faults/test_finite_fault_stack_wiring.py +++ b/tests/test_common/test_api/test_faults/test_finite_fault_stack_wiring.py @@ -9,7 +9,8 @@ _options_with_finite_fault_gradients, ) from gempy_engine.API.interp_single._multi_scalar_field_manager import _compute_independent_chunks -from gempy_engine.core.data import FiniteFault, InterpolationOptions +from gempy_engine.config import AvailableBackends +from gempy_engine.core.data import FiniteFault, InterpolationOptions, TaperType from gempy_engine.core.backend_tensor import BackendTensor from gempy_engine.core.data.exported_fields import ExportedFields from gempy_engine.core.data.kernel_classes.faults import FaultsData @@ -18,7 +19,7 @@ from gempy_engine.core.data.stacks_structure import StacksStructure -def test_modify_fault_values_applies_projected_taper(): +def test_modify_fault_values_reaches_zero_at_fault_tip(): finite_fault = FiniteFault(center=(0.0, 0.0, 0.0), strike_radius=1.0, dip_radius=1.0) fault_input = FaultsData.from_user_input(thickness=None, finite_fault=finite_fault) points = np.array([ @@ -47,6 +48,51 @@ def test_modify_fault_values_applies_projected_taper(): assert np.allclose(tapered_values, [[0.0, 2.0, 0.0]]) +@pytest.mark.parametrize("use_gpu", [False, True]) +def test_modify_fault_values_preserves_pytorch_backend(use_gpu): + torch = pytest.importorskip("torch") + if use_gpu and not torch.cuda.is_available(): + pytest.skip("CUDA is not available") + BackendTensor._change_backend(AvailableBackends.PYTORCH, use_gpu=use_gpu) + dtype = BackendTensor.dtype_obj + device = BackendTensor.device + + finite_fault = FiniteFault( + center=(0.0, 0.0, 0.0), + strike_radius=1.0, + dip_radius=1.0, + taper=TaperType.SPLINE, + ) + fault_input = FaultsData.from_user_input(thickness=None, finite_fault=finite_fault) + points = torch.tensor([ + [2.0, 0.0, 2.0], + [0.0, 0.0, 2.0], + [1.0, 0.0, 2.0], + ], dtype=dtype, device=device) + exported_fields = ExportedFields( + _scalar_field=torch.full((3,), 2.0, dtype=dtype, device=device), + _gx_field=torch.zeros(3, dtype=dtype, device=device), + _gy_field=torch.zeros(3, dtype=dtype, device=device), + _gz_field=torch.ones(3, dtype=dtype, device=device), + _grid_size=3, + _scalar_field_at_surface_points=torch.tensor([0.0], dtype=dtype, device=device), + ) + output = ScalarFieldOutput( + weights=None, + grid=None, + exported_fields=exported_fields, + stack_relation=StackRelationType.FAULT, + values_block=torch.tensor([[1.0, 3.0, 3.0]], dtype=dtype, device=device), + ) + + tapered_values = _modify_faults_values_output(fault_input, output, points) + + assert isinstance(tapered_values, torch.Tensor) + assert tapered_values.dtype == dtype + assert tapered_values.device == points.device + assert torch.allclose(tapered_values, torch.tensor([[0.0, 2.0, 0.0]], dtype=dtype, device=device)) + + def test_finite_fault_gradient_options_do_not_mutate_input(): options = InterpolationOptions.from_args(range=1.0, c_o=1.0, mesh_extraction=False) fault_input = FaultsData.from_user_input( diff --git a/tests/test_server/test_finite_fault_input.py b/tests/test_server/test_finite_fault_input.py new file mode 100644 index 00000000..07125601 --- /dev/null +++ b/tests/test_server/test_finite_fault_input.py @@ -0,0 +1,50 @@ +from gempy_engine.core.data.finite_fault import TaperType +from gempy_engine.core.data.input_data_descriptor import InputDataDescriptor +from gempy_engine.core.data.kernel_classes.server.input_parser import InputDataDescriptorSchema + + +def test_finite_fault_server_schema_converts_to_runtime_data(): + schema = InputDataDescriptorSchema.model_validate({ + "number_of_points_per_surface" : [3, 3], + "number_of_points_per_stack" : [3, 3], + "number_of_orientations_per_stack" : [1, 1], + "number_of_surfaces_per_stack" : [1, 1], + "masking_descriptor" : [3, 1], + "faults_relations" : [[False, True], [False, False]], + "faults_input_data" : [ + { + "finite_fault": { + "center" : [0.0, 0.0, 0.0], + "strike_radius": [2.0, 1.0], + "dip_radius" : 0.75, + "taper" : "quadratic", + "rotation_deg" : 15.0, + } + }, + None, + ], + }) + + descriptor = InputDataDescriptor.from_schema(schema) + fault_data = descriptor.stack_structure.faults_input_data[0] + + assert descriptor.stack_structure.faults_relations.tolist() == [[False, True], [False, False]] + assert fault_data.finite_fault.center == (0.0, 0.0, 0.0) + assert fault_data.finite_fault.strike_radius == (2.0, 1.0) + assert fault_data.finite_fault.taper is TaperType.QUADRATIC + + +def test_finite_fault_server_schema_serializes_as_json_values(): + schema = InputDataDescriptorSchema.model_validate({ + "number_of_points_per_surface" : [3], + "number_of_points_per_stack" : [3], + "number_of_orientations_per_stack" : [1], + "number_of_surfaces_per_stack" : [1], + "masking_descriptor" : [3], + "faults_input_data" : [{"finite_fault": {"center": [1.0, 2.0, 3.0]}}], + }) + + payload = schema.model_dump(mode="json") + + assert payload["faults_input_data"][0]["finite_fault"]["center"] == [1.0, 2.0, 3.0] + assert payload["faults_input_data"][0]["finite_fault"]["taper"] == "cubic"