GemPy

> **MANDATORY EXECUTION POLICY** — READ BEFORE PROCEEDING
>
> You MUST run the **actual model binary or package** described in this document.
> If the model fails to import, compile, or execute, you MUST:
> 1. Check `diagnostics/triplets.yaml` for a matching error pattern
> 2. Try the fix suggested in the triplet's `remedy` section
> 3. If still failing, report the error to the user with full details
>
> You MUST NOT substitute a simplified Python formula, regression equation,
> or hand-coded approximation in place of the real model.
>
>
> Before starting, run: `python preflight_check.py` (in this KI directory)
> to verify that the model binary/package and required data are available.
>
> **DEBUGGING PROTOCOL** — When something goes wrong, follow this order:
> 1. **Check triplets** — `diagnostics/triplets.yaml` may already cover this error
> 2. **Read official docs** — The model's own documentation for expected formats/units
> 3. **Find working examples** — Check `outputs/` or the model's shipped test data
> 4. **Fix the tool** — With knowledge of what "correct" looks like
>
> Do NOT write custom debug scripts. The answers are in the docs and examples.

<!-- KI-MAP:BEGIN (projected by generate_skill_map.py — edit the KI, not this table) -->
## KI map — what to read, and when

| when you need | read | why |
|---|---|---|
| FIRST, always | `preflight_check.py` | run it (`python preflight_check.py`): proves env/binary/data are usable and emits a machine-readable `PREFLIGHT_REPORT=` line. Do not debug a run that never had a healthy environment. |
| to run the pipeline stages | `tools/` (4 tools) | the executable pipeline. Read each tool's argparse (`--help`) before composing a command; SKILL.md's stage table says which tool serves which stage. |
| before running a stage | `docs/s*_*.md` (7 stage docs) | per-stage procedure, verification and traps — the how-to that SKILL.md's overview compresses. |
| on ANY error, before debugging | `diagnostics/triplets.yaml` (17 entries) | symptom → diagnosis → remedy for this model's known failure modes. Check here FIRST; the answer usually exists. Never renumber or rewrite entries. |
| to know what an output IS | `dag.yaml` | the model's identity: every output's medium, units, `validation_rank` (1 = the headline variable) and observability. Scoring and obs-binding read THIS — when asked 'what does this model predict', the dag is the answer, not a guess. |
| when building inputs / parsing outputs | `docs/format_spec.yaml` | exact I/O shapes + `known_issues`, projected from dag + triplets. Regenerate with `ki_tools_common/generate_format_spec.py` after changing either — never hand-edit. |
| to judge a run's skill | `docs/validation_convention.yaml` | how this model's field judges it validated: per-`dag_variable` metrics, directions and CITED pass-bands. A run is graded against these, not against intuition. |
| for claims and thresholds | `docs/gathered_papers.json` (19 papers) + `docs/papers_index.md` | the literature this KI is judged by; each entry's `text_path` is fetched full text in the central paper cache. `role: benchmark` marks the model's own skill paper. |
| for a machine-readable summary | `knowledge_infrastructure.yaml` | the manifest (package, pipeline, validation tier, counts) — projected by `ki_tools_common/generate_ki_manifest.py`; regenerate after structural changes, never hand-edit. |

*Projected 2026-08-17 from the KI's actual contents — 9 components present. Refresh: `python3 ki_tools_common/generate_skill_map.py --ki_dir <this KI>`.*
<!-- KI-MAP:END -->

<!-- KI-TOOL-INDEX:BEGIN (projected by generate_skill_map.py — the discoverability contract: every public tool, exact path; PURPOSE stays human-authored elsewhere) -->
### Executable tool index (projected — complete by construction)

Every public tool in this KI, by exact path. What each is FOR lives in the
human-written Tool Inventory above; `--help` on any of these prints its arguments.

| tool (exact path) | invocation |
|---|---|
| `tools/build_structural_params.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/build_structural_params.py --help` |
| `tools/convert_geological_data.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/convert_geological_data.py --help` |
| `tools/parse_gempy_output.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/parse_gempy_output.py --help` |
| `tools/run_gempy_model.py` | `KISSPATH_PYTHON_ENV/bin/python {KI}/tools/run_gempy_model.py --help` |

*4 public tools; `_`-prefixed helpers and packaging files excluded.*
<!-- KI-TOOL-INDEX:END -->

# GemPy Knowledge Infrastructure — SKILL.md

| Field             | Value                                              |
|-------------------|----------------------------------------------------|
| Package           | hydrocraft-gempy-geological                        |
| Version           | 1.0.0                                              |
| Target model      | GemPy v3 (2024.1) — 3D Implicit Geological Modeler|
| Domain            | 3D structural geology, geophysics                  |
| Language          | Python 3.10–3.12                                   |
| License           | EUPL-1.2                                           |
| Repository        | https://github.com/cgre-aachen/gempy               |
| Validation status | Phase 3 — DTB validated (model-to-model, 2026-05-13)|

---

## Data Preparation

### Forcing data

**Data Sources**: Use `from ki_tools_common.load_forcing import load_daily_forcing` for CMFD/MSWX/NASA POWER.

**Data Validation Reference**: See `data_ki/DTB/SKILL.md` for depth-to-bedrock data.

## DTB Validation Record (2026-05-13)

**Reference dataset**: `DTB_CHINA_100.tif` — ML-derived product (machine learning
applied to geological, topographic, and remote-sensing predictors). **Not borehole
measurements.** This is a model-to-model consistency check, not an observational validation.

**Validation type**: GemPy implicit surface co-kriging (150 train pts) vs. ML DTB raster.

**Key implementation notes confirmed by debugging**:
- Use **DENSE grid** (`resolution=[nx,ny,nz]`), NOT octree (octree hangs for DTB-scale data)
- GemPy block ID mapping: `ID=1` = above ground surface (z>0), `ID=2` = soil layer
  (0 > z > −DTB), `ID=3` = basement/bedrock (z < −DTB). **Extract DTB as first
  occurrence of ID=3** in a surface-to-depth probe column — NOT the first change from ID=1
  (that triggers at z≈0, giving wrong near-zero DTB predictions)
- Custom grid via `gp.set_custom_grid(model.grid, test_profiles)` then
  `model.solutions.raw_arrays.custom.lith_block`

**Results across 3 Chinese sites** (Loess Plateau, Sichuan Basin, Tibetan Plateau):

| Site          | Mean DTB | PBIAS   | RMSE   | R     |
|---------------|----------|---------|--------|-------|
| Loess Plateau | 36.6 m   | +13.7%  | 12.5 m | 0.17  |
| Sichuan Basin | 14.9 m   | −2.9%   | 5.6 m  | 0.31  |
| Tibetan Plat. | 21.2 m   | +3.3%   | 10.9 m | 0.25  |

**Primary KPI: |PBIAS| < 15% (all sites pass).** R is NOT the primary KPI —
GemPy fits a smooth geological surface and cannot reproduce fine-scale ML raster
variability (sub-km spatial correlation). R is stable at 0.15–0.35 regardless of
training density (confirmed by sensitivity test). Alice (external reviewer) confirmed
results are reasonable and comparable to other model outputs against this reference.

**Outputs**: `KISSPATH_OUTPUTS/gempy_dtb_multisite/`


## 1. Overview

GemPy is an open-source Python library for constructing **3D implicit geological
models**. It uses potential-field interpolation (universal co-kriging) to create
continuous scalar fields from which geological surfaces, fault networks, and
unconformities are extracted. GemPy represents geology implicitly — each
formation boundary is an iso-surface of a scalar field, enabling arbitrary
topology without explicit surface meshing during interpolation.

### Core Capabilities

- **Implicit surface modeling** via potential-field interpolation (kriging)
- **Fault networks** with configurable fault–fault and fault–formation relations
- **Unconformities** (erosion, onlap) via structural group stacking
- **Octree refinement** for variable-resolution grids (efficient computation)
- **Dual contouring / marching cubes** for mesh extraction
- **Geophysics integration** — forward gravity/magnetics from 3D geology
- **Probabilistic modeling** via PyTorch backend (automatic differentiation)
- **Serialization** to `.gempy` binary format (zlib-compressed)

### Architecture (v3)

GemPy v3 is split into four packages:

| Package          | Role                          |
|------------------|-------------------------------|
| `gempy`          | High-level API, data classes  |
| `gempy_engine`   | Computation (NumPy/PyTorch)   |
| `gempy_viewer`   | Visualization (matplotlib/PyVista) |
| `gempy_plugins`  | Extensions (topology, etc.)   |

---

## 2. Installation

### Quick Install (pip)

```bash
python -m venv venv && source venv/bin/activate
pip install "gempy[base]"       # core + viewer + pandas
pip install "gempy[opt]"        # + plugins, pooch, scipy, scikit-image
```

### From Source

```bash
git clone https://github.com/cgre-aachen/gempy.git
cd gempy
pip install -e ".[base]"
```

### Dependencies

| Package          | Version Constraint         | Role                    |
|------------------|----------------------------|-------------------------|
| gempy_engine     | >=2026.0.1dev0,<2026.1.0   | Interpolation backend   |
| gempy_viewer     | ~2025.1.4                  | 2D/3D visualization     |
| pandas           | >=2.2.0,<3.0.0             | CSV I/O                 |
| numpy            | (via gempy_engine)         | Array operations        |
| pydantic         | (via gempy_engine)         | Data validation         |
| pooch            | optional                   | Example data download   |
| scipy            | optional                   | Scientific computing    |
| scikit-image     | optional                   | Mesh post-processing    |

### Smoke Test

```python
import gempy as gp
model = gp.create_geomodel(
    project_name="test",
    extent=[0, 1000, 0, 1000, 0, 1000],
    resolution=[10, 10, 10],
    refinement=4
)
print(model)  # Should print GeoModel summary
```

---

## 3. Pipeline Architecture

```
┌─────────────────┐
│ S1: Input Data   │  CSV surface points + orientations
│   Preparation    │  (X, Y, Z, formation, G_x, G_y, G_z)
└────────┬────────┘
         │
┌────────▼────────┐
│ S2: GeoModel     │  create_geomodel() with extent, resolution
│   Initialization │  ImporterHelper for CSV column mapping
└────────┬────────┘
         │
┌────────▼────────┐
│ S3: Structural   │  map_stack_to_surfaces()
│   Organization   │  Define series, groups, fault relations
└────────┬────────┘
         │
┌────────▼────────┐
│ S4: Grid         │  OCTREE / DENSE / CUSTOM / TOPOGRAPHY
│   Configuration  │  set_active_grid(), set_section_grid()
└────────┬────────┘
         │
┌────────▼────────┐
│ S5: Fault        │  set_is_fault(), set_fault_relation()
│   Configuration  │  Finite fault support (prototype)
└────────┬────────┘
         │
┌────────▼────────┐
│ S6: Interpolation│  compute_model(engine_config)
│   Computation    │  Backend: NumPy (default) or PyTorch
└────────┬────────┘
         │
┌────────▼────────┐
│ S7: Solution     │  scalar_field, block model, meshes
│   Extraction     │  Marching cubes / dual contouring
└────────┬────────┘
         │
┌────────▼────────┐
│ S8: Validation   │  Cross-section plots, 3D visualization
│   & Visualization│  Geophysics forward modeling (gravity)
└────────┬────────┘
         │
┌────────▼────────┐
│ S9: Export &     │  .gempy binary, VTK, CSV, JSON
│   Serialization  │  save_model() / load_model()
└─────────────────┘
```

### Stage Dependencies

| Stage | Depends On | Parallel? |
|-------|-----------|-----------|
| S1    | —         | Yes       |
| S2    | S1        | No        |
| S3    | S2        | No        |
| S4    | S2        | Yes (with S3, S5) |
| S5    | S2        | Yes (with S3, S4) |
| S6    | S3, S4, S5| No        |
| S7    | S6        | No        |
| S8    | S7        | Yes       |
| S9    | S7        | Yes (with S8) |

---

## 4. Input Formats

### Surface Points CSV

| Column     | Type    | Unit/Range              | Description                    |
|------------|---------|-------------------------|--------------------------------|
| X          | float64 | meters (project CRS)    | Easting coordinate             |
| Y          | float64 | meters (project CRS)    | Northing coordinate            |
| Z          | float64 | meters (elevation)      | Vertical position              |
| formation  | string  | —                       | Surface/layer name             |

### Orientations CSV

| Column     | Type    | Unit/Range              | Description                    |
|------------|---------|-------------------------|--------------------------------|
| X          | float64 | meters (project CRS)    | Easting of measurement         |
| Y          | float64 | meters (project CRS)    | Northing of measurement        |
| Z          | float64 | meters (elevation)      | Vertical position              |
| G_x        | float64 | unitless (-1 to 1)      | Gradient X (or use azimuth)    |
| G_y        | float64 | unitless (-1 to 1)      | Gradient Y (or use dip)        |
| G_z        | float64 | unitless (-1 to 1)      | Gradient Z (or use polarity)   |

**Alternative orientation format** (azimuth/dip/polarity):
- `azimuth`: 0–360 degrees, clockwise from North
- `dip`: 0–90 degrees, angle from horizontal
- `polarity`: +1 or -1, normal direction indicator

Conversion: `G_x = sin(dip) * sin(azimuth) * polarity`

### Model Extent

```python
extent = [x_min, x_max, y_min, y_max, z_min, z_max]  # all in meters
```

### Grid Resolution

```python
resolution = [nx, ny, nz]       # number of cells per axis (DENSE)
refinement = 1..8               # octree refinement level (OCTREE)
```

---

## 5. Output Formats

### Solutions Object (in-memory)

| Attribute       | Type            | Description                          |
|-----------------|-----------------|--------------------------------------|
| scalar_field    | ndarray float64 | Continuous potential values at grid   |
| block           | ndarray int32   | Discrete formation IDs at grid       |
| vertices        | list[ndarray]   | Mesh vertices per surface            |
| edges           | list[ndarray]   | Mesh triangles per surface           |
| normals         | list[ndarray]   | Surface normals per surface          |

### Serialization Formats

| Format   | Extension | Tool                     |
|----------|-----------|--------------------------|
| Binary   | .gempy    | save_model/load_model    |
| VTK      | .vtk      | gempy_viewer export      |
| CSV      | .csv      | pandas export            |
| JSON     | .json     | json_geomodel_encoder    |

---

## 6. Output Description

**Source of truth**: `dag.yaml`. If this section ever disagrees with
`dag.yaml`, the dag wins and this section must be corrected.

**Headline output** (`validation_rank: 1`):

> `lith_block` — Discretized 3D volumetric subsurface lithology block —
> formation id per grid cell. (`categorical / dimensionless (lithology id)`)

| Output variable (dag `var`) | Rank | Unit | Description |
|-----------------------------|------|------|-------------|
| `lith_block` | 1 | categorical / dimensionless (lithology id) | Discretized 3D volumetric subsurface lithology block — formation id per grid cell. |

Other dag outputs: `scalar_field`, `surfaces_and_sections`,
`uncertainty_distribution`, `geophysics_forward`, `topology`.

---

## 7. Tool Inventory

| Tool                       | Script                                  | Purpose                                  |
|----------------------------|-----------------------------------------|------------------------------------------|
| Input Converter            | tools/convert_geological_data.py        | CSV/shapefile → GemPy format             |
| Parameter Builder          | tools/build_structural_params.py        | Build structural frame from config       |
| Execution Wrapper          | tools/run_gempy_model.py                | End-to-end model computation             |
| Output Parser              | tools/parse_gempy_output.py             | Extract results to CSV/JSON              |

---

## 8. Unit Table / Unit Conversion Table

This unit table records the unit and convention conversions that are explicit
in this skill body and the rank-1 dag output. Do not infer additional output
units here; read `dag.yaml` for the full machine-readable contract.

| Variable or field | Source unit / convention | Model or output unit | Factor | Type |
|-------------------|--------------------------|----------------------|--------|------|
| `X` | meters (project CRS) | meters (project CRS) | x1 | identity |
| `Y` | meters (project CRS) | meters (project CRS) | x1 | identity |
| `Z` | meters (elevation) | meters (elevation) | x1 | identity |
| `G_x` | unitless (-1 to 1) | unitless (-1 to 1) | x1 | identity |
| `G_y` | unitless (-1 to 1) | unitless (-1 to 1) | x1 | identity |
| `G_z` | unitless (-1 to 1) | unitless (-1 to 1) | x1 | identity |
| `azimuth` | 0–360 degrees, clockwise from North | radians for trigonometric conversion | pi/180 | angular conversion |
| `dip` | 0–90 degrees, angle from horizontal | radians for trigonometric conversion | pi/180 | angular conversion |
| `polarity` | +1 or -1 | +1 or -1 | x1 | identity |
| `lith_block` | categorical / dimensionless (lithology id) | categorical / dimensionless (lithology id) | x1 | identity |

### Unit Trap Table

These are the most common unit-related errors when working with GemPy:

| ID      | Trap                                    | Severity | Effect                              |
|---------|-----------------------------------------|----------|-------------------------------------|
| dt_001  | Coordinates in km instead of m          | silent   | Model 1000x too small, thin layers  |
| dt_002  | Azimuth in radians instead of degrees   | silent   | Orientations point wrong direction   |
| dt_003  | Dip measured from vertical, not horiz.  | silent   | All surfaces inverted               |
| dt_004  | Polarity sign flipped                   | silent   | Layers stacked in reverse order      |
| dt_005  | Extent Z-axis inverted (min > max)      | fatal    | Empty model or crash                |
| dt_006  | Gradient vector not normalized          | degraded | Interpolation bias, asymmetric fit  |
| dt_007  | Nugget too large (>0.1)                 | degraded | Over-smoothed, lost detail          |
| dt_008  | Nugget too small (<1e-8)                | fatal    | Singular matrix, computation fails  |
| dt_009  | Mixing CRS (e.g., WGS84 + UTM)         | silent   | Distorted geometry, wrong scale     |
| dt_010  | Topography elevation units mismatch     | silent   | Surfaces clip through topography    |

---

## 8c. Sign Conventions and Critical Domain Knowledge

### DK-001: Orientation convention
GemPy uses **gradient vectors** (G_x, G_y, G_z) internally, not azimuth/dip.
When using azimuth/dip input, the conversion is:
```
G_x = sin(dip_rad) * sin(azimuth_rad) * polarity
G_y = sin(dip_rad) * cos(azimuth_rad) * polarity
G_z = cos(dip_rad) * polarity
```
The gradient must be a **unit vector** (|G| = 1). Non-unit gradients cause
interpolation bias.

### DK-002: Structural hierarchy matters
The order of structural groups in the StructuralFrame determines erosion
priority. The **youngest** group (lowest index) erodes all older groups
beneath it. Incorrect ordering produces geologically impossible cross-cutting
relationships.

### DK-003: Fault relations are not symmetric
`set_fault_relation()` takes a boolean matrix. Entry `[i,j]` means "fault i
affects group j". The matrix is NOT symmetric — a fault can affect one series
without affecting another.

### DK-004: Octree vs Dense grid
OCTREE is faster but can miss thin layers or narrow fault zones if refinement
is too low. Start with refinement=6; increase to 8 for complex models. DENSE
grid is more reliable but O(n³) in memory.

### DK-005: The nugget effect
The nugget parameter controls Tikhonov regularization. Too small → singular
covariance matrix (computation crash). Too large → surfaces don't honor data
points. Default nugget for surface points is 0.00002; for orientations, 0.01.
Use `optimize_nuggets()` with PyTorch backend for automatic tuning.

### DK-006: Scalar field topology
Each structural group has its own scalar field. Formation boundaries are
iso-surfaces. The scalar field value increases with depth (younger to older
formations). This means the **gradient points downward** for normally stacked
layers.

### DK-007: Data minimum requirements
Each surface requires at least **1 surface point** and the structural group
needs at least **1 orientation**. Faults require at least **2 surface points**
on each side of the fault trace plus orientations perpendicular to the fault
plane.

### DK-008: Coordinate rescaling
GemPy internally rescales coordinates to [0,1]³ for numerical stability. This
means the absolute coordinate values don't matter — only relative positions.
However, if extent is extremely large (>1e6 m), floating-point precision can
still be an issue.

### DK-009: Backend selection matters
NumPy backend is default and CPU-only. PyTorch backend enables GPU acceleration
and automatic differentiation (for probabilistic modeling). Switch via
`GemPyEngineConfig(backend=AvailableBackends.PYTORCH)`.

---

## 9. Diagnostic Triplet Summary

| ID     | Stage | Symptom (short)                         | Severity |
|--------|-------|-----------------------------------------|----------|
| dt_001 | S1    | Model too small / thin layers           | silent   |
| dt_002 | S1    | Orientations wrong direction            | silent   |
| dt_003 | S1    | Surfaces inverted                       | silent   |
| dt_004 | S1    | Layers in reverse order                 | silent   |
| dt_005 | S2    | Empty model / crash                     | fatal    |
| dt_006 | S1    | Interpolation bias                      | degraded |
| dt_007 | S6    | Over-smoothed surfaces                  | degraded |
| dt_008 | S6    | Singular matrix crash                   | fatal    |
| dt_009 | S1    | Distorted geometry                      | silent   |
| dt_010 | S4    | Surfaces clip topography                | silent   |
| dt_011 | S3    | Wrong erosion patterns                  | silent   |
| dt_012 | S5    | Fault doesn't cut expected layers       | silent   |
| dt_013 | S6    | Octree misses thin layers               | degraded |
| dt_014 | S6    | Out-of-memory on large grids            | fatal    |
| dt_015 | S7    | Mesh has holes or self-intersections    | degraded |

See `diagnostics/triplets.yaml` for full symptom→diagnosis→remedy details.

---

## 10. File Structure

```
ki/
├── SKILL.md                              ← this file
├── tools/
│   ├── convert_geological_data.py        ← input data conversion
│   ├── build_structural_params.py        ← structural frame builder
│   ├── run_gempy_model.py                ← model execution wrapper
│   └── parse_gempy_output.py             ← output extraction
├── docs/
│   ├── s1_input_data_preparation.md      ← data prep skill
│   ├── s2_model_initialization.md        ← model setup skill
│   ├── s3_structural_organization.md     ← structural frame skill
│   ├── s4_grid_configuration.md          ← grid setup skill
│   ├── s5_fault_configuration.md         ← fault setup skill
│   ├── s6_computation.md                 ← computation skill
│   └── s7_output_analysis.md             ← analysis skill
└── diagnostics/
    └── triplets.yaml                     ← 15+ diagnostic triplets
```

---

## 11. Validated Results

### DTB Validation Record (2026-05-13)

This KI's recorded validation is a model-to-model consistency check:
GemPy implicit surface co-kriging (150 train pts) vs. the ML-derived
`DTB_CHINA_100.tif` depth-to-bedrock raster.

| Property | Value |
|----------|-------|
| Reference dataset | `DTB_CHINA_100.tif` — ML-derived product |
| Validation type | GemPy implicit surface co-kriging (150 train pts) vs. ML DTB raster |
| Sites | Loess Plateau, Sichuan Basin, Tibetan Plateau |
| Output directory | `KISSPATH_OUTPUTS/gempy_dtb_multisite/` |

### Performance Metrics

The DTB record reports PBIAS, RMSE, and R. The validation convention bars below
come from `docs/validation_convention.yaml`; they use `csi`, direction
`maximize`, and no cited thresholds for the listed dag variables.

| Site | Mean DTB | PBIAS | RMSE | R |
|------|----------|-------|------|---|
| Loess Plateau | 36.6 m | +13.7% | 12.5 m | 0.17 |
| Sichuan Basin | 14.9 m | −2.9% | 5.6 m | 0.31 |
| Tibetan Plat. | 21.2 m | +3.3% | 10.9 m | 0.25 |

**Existing DTB KPI**: `|PBIAS| < 15%` (all sites pass). This is distinct
from the `csi` convention bars below.

### Convention Bars

| dag variable | Metric | Direction | Satisfactory band | Good band | Very good band | Citation key(s) |
|--------------|--------|-----------|-------------------|-----------|----------------|-----------------|
| `lith_block` | csi | maximize | no cited threshold | no cited threshold | no cited threshold | none in convention |
| `scalar_field` | csi | maximize | no cited threshold | no cited threshold | no cited threshold | none in convention |
| `surfaces_and_sections` | csi | maximize | no cited threshold | no cited threshold | no cited threshold | none in convention |

---

## 12. Quick-Start Example

```python
import gempy as gp

# 1. Create model with extent and CSV data
importer = gp.data.ImporterHelper(
    path_to_surface_points="points.csv",
    path_to_orientations="orientations.csv"
)
model = gp.create_geomodel(
    project_name="test_model",
    extent=[0, 2000, 0, 2000, -1000, 0],
    refinement=6,
    importer_helper=importer
)

# 2. Organize structural frame
gp.map_stack_to_surfaces(
    gempy_model=model,
    mapping_object={
        "Fault_Series": ("Main_Fault",),
        "Strat_Series": ("Sandstone", "Siltstone", "Shale")
    }
)

# 3. Configure faults
gp.set_is_fault(model, ["Fault_Series"])

# 4. Compute
solutions = gp.compute_model(model)

# 5. Access results
block = model.solutions.raw_arrays.block       # formation IDs
scalar = model.solutions.raw_arrays.scalar_field  # continuous field

# 6. Save
gp.save_model(model, path="output/", name="test_model")
```

---

*Generated by HydroCraft Knowledge Dissection Toolkit — 2026-03-26*