use crate::error::StaticError;
use crate::gridder::{
layer_grid_stack, Conformity, Control, ExtrapolationPolicy, KernelSurface, SolveOpts,
StreamingLayering, Surface, ZoneLayerSpec,
};
use crate::model::model::{Georef, StaticModel};
use crate::model::pipeline::{PropertyPipeline, PropertyReport};
use crate::model::population::{override_zone_priors, populate, PetroSample};
use crate::model::provenance::{
BuildWarning, HorizonTieResidual, InterfaceRepair, PopulationMode, Provenance, StackProvenance,
WellTieRecord, ZoneProvenance,
};
use crate::model::spec::{BuildSpec, TieMethod, TieSettings};
use crate::model::trend::TrendSurface;
use crate::model::zones::ZoneTable;
use crate::spill::{decide_mode, BuildMode, MemoryBudget, SpillNotice};
use crate::volumetrics::{validate_positive, ConstantPriors, NTG, PORO, SW};
use crate::wireframe::{
Contact, ContactKind, GriddedDepth, Hardness, Horizon, HorizonRole, Wireframe,
};
use petektools::{grid_min_curvature_conditioned, Conditioning, Lattice, MinCurvatureOperator};
use rayon::prelude::*;
use serde::{Deserialize, Serialize};
use std::path::PathBuf;
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct BuildOpts {
pub area_m2: f64,
pub gross_height_m: f64,
pub nk: usize,
pub conformity: Conformity,
pub solve_opts: SolveOpts,
pub priors: ConstantPriors,
}
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct Pick {
pub ip: usize,
pub jp: usize,
pub depth_m: f64,
}
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct WorldPoint {
pub x: f64,
pub y: f64,
pub depth_m: f64,
}
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub enum HorizonSource {
Scatter(Vec<WorldPoint>),
Mapped(GriddedDepth),
TopsOnly(Vec<Pick>),
}
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct WellTie {
pub id: String,
pub x: f64,
pub y: f64,
pub ip: usize,
pub jp: usize,
pub tops: Vec<(String, f64)>,
}
impl WellTie {
#[must_use]
pub fn new(id: impl Into<String>, x: f64, y: f64, ip: usize, jp: usize) -> Self {
Self {
id: id.into(),
x,
y,
ip,
jp,
tops: Vec::new(),
}
}
#[must_use]
pub fn with_top(mut self, horizon: impl Into<String>, depth_m: f64) -> Self {
let horizon = horizon.into();
self.tops.retain(|(h, _)| h != &horizon);
self.tops.push((horizon, depth_m));
self
}
}
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct StackHorizon {
pub name: String,
pub source: HorizonSource,
}
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct StackZone {
pub name: String,
pub color: Option<String>,
pub conformity: Conformity,
pub nk: usize,
pub contacts: Vec<Contact>,
}
impl StackZone {
#[must_use]
pub fn new(
name: impl Into<String>,
conformity: Conformity,
nk: usize,
contacts: Vec<Contact>,
) -> Self {
Self {
name: name.into(),
color: None,
conformity,
nk,
contacts,
}
}
#[must_use]
pub fn with_color(mut self, color: impl Into<String>) -> Self {
self.color = Some(color.into());
self
}
}
#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
pub struct HorizonStack {
pub horizons: Vec<StackHorizon>,
pub zone_layers: Vec<StackZone>,
}
#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
pub struct StackFrame {
pub ni: usize,
pub nj: usize,
pub georef: Georef,
}
#[derive(Debug, Clone)]
pub struct StaticModelBuilder {
opts: BuildOpts,
ni: usize,
nj: usize,
controls: Vec<Control>,
base_controls: Option<Vec<Control>>,
logs: Option<Vec<PetroSample>>,
trend: Option<TrendSurface>,
framework: Wireframe,
default_inputs_ref: &'static str,
spec: BuildSpec,
warnings: Vec<BuildWarning>,
properties: Vec<PropertyPipeline>,
stack: Option<HorizonStack>,
zone_properties: Vec<(String, PropertyPipeline)>,
zone_priors: Vec<(String, ConstantPriors)>,
budget: MemoryBudget,
spill_dir: Option<PathBuf>,
spill_persist: bool,
}
impl StaticModelBuilder {
pub fn flat(
ni: usize,
nj: usize,
top_depth_m: f64,
contact_depth_m: f64,
opts: BuildOpts,
) -> Result<Self, StaticError> {
if ni == 0 || nj == 0 || opts.nk == 0 {
return Err(StaticError::Grid("dimensions must be non-zero".into()));
}
validate_positive("area_m2", opts.area_m2)?;
validate_positive("gross_height_m", opts.gross_height_m)?;
if !contact_depth_m.is_finite() {
return Err(StaticError::InvalidInput(format!(
"contact depth must be finite, got {contact_depth_m}"
)));
}
let controls: Vec<Control> = [(0, 0), (ni, 0), (0, nj), (ni, nj)]
.iter()
.map(|&(ip, jp)| Control {
ip,
jp,
z: top_depth_m,
})
.collect();
let framework = synth_framework(ni + 1, nj + 1, top_depth_m, contact_depth_m);
Ok(Self {
opts,
ni,
nj,
controls,
base_controls: None,
logs: None,
trend: None,
framework,
default_inputs_ref: "flat-box",
spec: BuildSpec::default(),
warnings: Vec::new(),
properties: Vec::new(),
stack: None,
zone_properties: Vec::new(),
zone_priors: Vec::new(),
budget: MemoryBudget::default(),
spill_dir: None,
spill_persist: false,
})
}
pub fn from_wireframe(wf: &Wireframe, opts: BuildOpts) -> Result<Self, StaticError> {
let s = wireframe_structure(wf, &opts)?;
Ok(Self {
opts,
ni: s.ni,
nj: s.nj,
controls: s.top_controls,
base_controls: s.base_controls,
logs: None,
trend: None,
framework: wf.clone(),
default_inputs_ref: "wireframe",
spec: BuildSpec::default(),
warnings: s.warnings,
properties: Vec::new(),
stack: None,
zone_properties: Vec::new(),
zone_priors: Vec::new(),
budget: MemoryBudget::default(),
spill_dir: None,
spill_persist: false,
})
}
pub fn from_horizon_stack(stack: HorizonStack, opts: BuildOpts) -> Result<Self, StaticError> {
let (ni, nj) = validate_stack(&stack)?;
Ok(Self {
opts,
ni,
nj,
controls: Vec::new(),
base_controls: None,
logs: None,
trend: None,
framework: Wireframe::from_boundary(crate::wireframe::Boundary {
ring: vec![[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]],
hardness: Hardness::Interpolated,
}),
default_inputs_ref: "horizon-stack",
spec: BuildSpec::default(),
warnings: Vec::new(),
properties: Vec::new(),
stack: Some(stack),
zone_properties: Vec::new(),
zone_priors: Vec::new(),
budget: MemoryBudget::default(),
spill_dir: None,
spill_persist: false,
})
}
pub fn from_scatter_stack(
mut stack: HorizonStack,
opts: BuildOpts,
frame: StackFrame,
) -> Result<Self, StaticError> {
condition_scatter(&mut stack, &frame)?;
let g = frame.georef;
Ok(Self::from_horizon_stack(stack, opts)?.with_georef(
g.origin_x,
g.origin_y,
g.spacing_x,
g.spacing_y,
))
}
pub fn condition_scatter_stack(
mut stack: HorizonStack,
frame: &StackFrame,
) -> Result<HorizonStack, StaticError> {
condition_scatter(&mut stack, frame)?;
Ok(stack)
}
#[must_use]
pub fn with_clamp_base_to_top(mut self, clamp: bool) -> Self {
self.spec.clamp_base_to_top = clamp;
self
}
#[must_use]
pub fn with_min_thickness_m(mut self, min_thickness_m: f64) -> Self {
self.spec.min_thickness_m = Some(min_thickness_m);
self
}
#[must_use]
pub fn with_collapse_below_m(mut self, collapse_below_m: f64) -> Self {
self.spec.collapse_below_m = Some(collapse_below_m);
self
}
#[must_use]
pub fn with_extrapolation(mut self, policy: ExtrapolationPolicy) -> Self {
self.spec.extrapolation = policy;
self
}
#[must_use]
pub fn with_sw_gas(mut self, sw_gas: f64) -> Self {
self.spec.sw_gas = Some(sw_gas);
self
}
#[must_use]
pub fn with_logs(mut self, samples: Vec<PetroSample>) -> Self {
self.logs = if samples.is_empty() {
None
} else {
Some(crate::model::population::sort_by_tvd(samples))
};
self
}
#[must_use]
pub fn with_areal_trend(mut self, trend: TrendSurface) -> Self {
self.trend = Some(trend);
self
}
#[must_use]
pub fn with_inputs_ref(mut self, inputs_ref: impl Into<String>) -> Self {
self.spec.inputs_ref = Some(inputs_ref.into());
self
}
#[must_use]
pub fn with_spec(mut self, spec: BuildSpec) -> Self {
self.spec = spec;
self
}
fn inputs_ref_string(&self) -> String {
self.spec
.inputs_ref
.clone()
.unwrap_or_else(|| self.default_inputs_ref.to_string())
}
#[must_use]
pub fn with_memory_budget(mut self, budget: MemoryBudget) -> Self {
self.budget = budget;
self
}
#[must_use]
pub fn with_spill_dir(mut self, dir: impl Into<PathBuf>) -> Self {
self.spill_dir = Some(dir.into());
self
}
#[must_use]
pub fn with_spill_persist(mut self, persist: bool) -> Self {
self.spill_persist = persist;
self
}
#[must_use]
pub fn with_georef(
mut self,
origin_x: f64,
origin_y: f64,
spacing_x: f64,
spacing_y: f64,
) -> Self {
self.spec.georef = Georef::new(origin_x, origin_y, spacing_x, spacing_y);
self
}
#[must_use]
pub fn with_boundary(mut self, ring: Vec<[f64; 2]>) -> Self {
self.spec.boundary = if ring.is_empty() { None } else { Some(ring) };
self
}
#[must_use]
pub fn with_property(mut self, pipeline: PropertyPipeline) -> Self {
self.properties.push(pipeline);
self
}
#[must_use]
pub fn with_zone_property(
mut self,
zone_name: impl Into<String>,
pipeline: PropertyPipeline,
) -> Self {
self.zone_properties.push((zone_name.into(), pipeline));
self
}
#[must_use]
pub fn with_zone_priors(
mut self,
zone_name: impl Into<String>,
priors: ConstantPriors,
) -> Self {
self.zone_priors.push((zone_name.into(), priors));
self
}
#[must_use]
pub fn with_tie_settings(mut self, ties: TieSettings) -> Self {
self.spec.ties = ties;
self
}
#[must_use]
pub fn with_well_ties(mut self, ties: Vec<WellTie>) -> Self {
self.spec.well_ties = ties;
self
}
#[must_use]
pub fn with_sugar_cube(mut self, sugar_cube: bool) -> Self {
self.spec.sugar_cube = sugar_cube;
self
}
pub fn add_top_control(&mut self, ip: usize, jp: usize, depth_m: f64) {
self.controls.push(Control { ip, jp, z: depth_m });
}
#[must_use]
pub fn control_count(&self) -> usize {
self.controls.len()
}
pub fn build(&self) -> Result<StaticModel, StaticError> {
if self.stack.is_some() {
return self.build_stack();
}
let nx = self.ni + 1;
let ny = self.nj + 1;
let top_k = flat_corner_surface(nx, ny, &self.controls)
.unwrap_or_else(|| warm_surface(nx, ny, &self.controls))?;
let top = top_k.surface();
let base = match &self.base_controls {
Some(bc) => warm_surface(nx, ny, bc)?.surface().clone(),
None => top.offset_by(self.opts.gross_height_m),
};
let mut warnings = self.warnings.clone();
let base = match self.spec.min_thickness_m {
Some(min_t) => {
let (repaired, columns, worst_m) = base.repair_min_thickness(top, min_t)?;
if columns > 0 {
warnings.push(BuildWarning::ThinColumnsRepaired { columns, worst_m });
}
repaired
}
None => base.guard_below(top, self.spec.clamp_base_to_top)?,
};
let (dx, dy) = spacing(self.opts.area_m2, self.ni, self.nj);
if self.streaming_eligible() {
let streaming = StreamingLayering::prepare(
&[top, &base],
dx,
dy,
&[ZoneLayerSpec {
conformity: self.opts.conformity,
requested_nk: self.opts.nk,
}],
)?;
let (mode, estimate) = decide_mode(streaming.dims(), 3, self.budget);
if mode == BuildMode::Spilled {
return self.build_spilled_streaming(&streaming, estimate, warnings);
}
}
let layered = layer_grid_stack(
&[top, &base],
dx,
dy,
&[ZoneLayerSpec {
conformity: self.opts.conformity,
requested_nk: self.opts.nk,
}],
self.spec.collapse_below_m,
)?;
let nk = layered.nk;
if layered.truncated_cells > 0 {
warnings.push(BuildWarning::LayersTruncated {
cells: layered.truncated_cells,
});
}
if layered.collapsed_cells > 0 {
warnings.push(BuildWarning::CellsCollapsed {
cells: layered.collapsed_cells,
});
}
if layered.nk_capped {
warnings.push(BuildWarning::LayerCountCapped { nk });
}
let mut grid = layered.grid;
populate(
&mut grid,
self.opts.priors,
self.logs.as_deref(),
self.trend.as_ref(),
)?;
let mut property_reports = Vec::with_capacity(self.properties.len());
for pipe in &self.properties {
property_reports.push(pipe.apply_with_georef(&mut grid, self.spec.georef)?);
}
let zone_kranges = vec![("RESERVOIR".to_string(), 0..nk)];
self.apply_zone_population(&mut grid, &zone_kranges, &mut property_reports)?;
let population = self.population_mode();
self.maybe_spill(
StaticModel::new(
self.framework.clone(),
grid,
ZoneTable::single(nk),
Provenance {
inputs_ref: self.inputs_ref_string(),
solve_opts: self.opts.solve_opts,
conformity: self.opts.conformity,
nk,
population,
realization: None,
warnings,
property_reports,
stack: None,
well_ties: Vec::new(),
sugar_cube: self.spec.sugar_cube,
},
self.spec.sw_gas,
)
.with_georef_opt(self.spec.georef),
)
}
fn maybe_spill(&self, model: StaticModel) -> Result<StaticModel, StaticError> {
let dims = model.grid().dims();
let n_cubes = model.property_names().len();
let (mode, estimate) = decide_mode(dims, n_cubes, self.budget);
if mode == BuildMode::InCore {
return Ok(model);
}
let dir = self.spill_dir.clone().unwrap_or_else(std::env::temp_dir);
let backing = crate::spill::spill_grid(model.grid(), &dir, !self.spill_persist)?;
SpillNotice {
cells: dims.cell_count(),
budget_bytes: self.budget.limit_bytes(),
estimate_bytes: estimate,
store_path: backing.store_path().to_path_buf(),
}
.warn();
Ok(model.into_spilled(backing))
}
fn streaming_eligible(&self) -> bool {
self.spec.collapse_below_m.is_none()
&& self.logs.is_none()
&& self.trend.is_none()
&& self.properties.is_empty()
&& self.zone_properties.is_empty()
}
fn build_spilled_streaming(
&self,
streaming: &StreamingLayering,
estimate: u64,
mut warnings: Vec<BuildWarning>,
) -> Result<StaticModel, StaticError> {
let dims = streaming.dims();
let nk = dims.nk;
let mut cube_names = vec![NTG.to_string(), PORO.to_string(), SW.to_string()];
cube_names.sort();
let priors = self.opts.priors;
let const_for = |name: &str| -> f32 {
(if name == PORO {
priors.porosity
} else if name == NTG {
priors.net_to_gross
} else {
priors.water_saturation
}) as f32
};
let mut coord = Vec::new();
streaming.fill_coord(&mut coord);
let plane_len = dims.pillar_count();
let mut plane_top = vec![0.0f64; plane_len];
let mut plane_bot = vec![0.0f64; plane_len];
let mut truncated = 0usize;
let dir = self.spill_dir.clone().unwrap_or_else(std::env::temp_dir);
let path = crate::spill::unique_spill_path(&dir);
let backing = crate::spill::spill_streaming(
&path,
dims,
&coord,
&cube_names,
!self.spill_persist,
|k, out| {
truncated += streaming.fill_zcorn_slab(k, &mut plane_top, &mut plane_bot, out);
Ok(())
},
|name, _k, out| {
out.fill(const_for(name));
Ok(())
},
)?;
if truncated > 0 {
warnings.push(BuildWarning::LayersTruncated { cells: truncated });
}
if streaming.nk_capped() {
warnings.push(BuildWarning::LayerCountCapped { nk });
}
SpillNotice {
cells: dims.cell_count(),
budget_bytes: self.budget.limit_bytes(),
estimate_bytes: estimate,
store_path: backing.store_path().to_path_buf(),
}
.warn();
Ok(StaticModel::spilled(
self.framework.clone(),
ZoneTable::single(nk),
Provenance {
inputs_ref: self.inputs_ref_string(),
solve_opts: self.opts.solve_opts,
conformity: self.opts.conformity,
nk,
population: self.population_mode(),
realization: None,
warnings,
property_reports: Vec::new(),
stack: None,
well_ties: Vec::new(),
sugar_cube: self.spec.sugar_cube,
},
self.spec.sw_gas,
self.spec.georef,
backing,
))
}
fn population_mode(&self) -> PopulationMode {
if self.logs.is_some() || !self.properties.is_empty() || !self.zone_properties.is_empty() {
PopulationMode::Logs
} else {
PopulationMode::Priors
}
}
fn apply_zone_population(
&self,
grid: &mut crate::grid::Grid,
zone_kranges: &[(String, core::ops::Range<usize>)],
reports: &mut Vec<PropertyReport>,
) -> Result<(), StaticError> {
let find = |name: &str| -> Result<core::ops::Range<usize>, StaticError> {
zone_kranges
.iter()
.find(|(n, _)| n == name)
.map(|(_, r)| r.clone())
.ok_or_else(|| {
StaticError::InvalidInput(format!(
"per-zone population: zone '{name}' is not in the model's zone table"
))
})
};
for (name, priors) in &self.zone_priors {
override_zone_priors(grid, *priors, find(name)?)?;
}
for (name, pipe) in &self.zone_properties {
let report = pipe
.apply_in_zone_with_georef(grid, find(name)?, self.spec.georef)
.map_err(|e| StaticError::InvalidInput(format!("zone '{name}': {e}")))?;
reports.push(report);
}
Ok(())
}
fn build_stack(&self) -> Result<StaticModel, StaticError> {
let stack = self.stack.as_ref().expect("build_stack requires a stack");
let nx = self.ni + 1;
let ny = self.nj + 1;
let mut surfaces = resolve_stack_surfaces(stack, nx, ny, self.spec.extrapolation)?;
let is_derived: Vec<bool> = stack
.horizons
.iter()
.map(|h| matches!(h.source, HorizonSource::TopsOnly(_)))
.collect();
let well_ties = self.apply_well_ties(stack, &mut surfaces, nx, ny)?;
let mut interface_repairs = Vec::new();
for i in 0..surfaces.len() - 1 {
let (upper, lower) = repair_interface(
&surfaces[i],
&surfaces[i + 1],
self.spec.min_thickness_m,
self.spec.clamp_base_to_top,
is_derived[i],
is_derived[i + 1],
i,
&mut interface_repairs,
)?;
surfaces[i] = upper;
surfaces[i + 1] = lower;
}
let specs: Vec<ZoneLayerSpec> = stack
.zone_layers
.iter()
.map(|z| ZoneLayerSpec {
conformity: z.conformity,
requested_nk: z.nk,
})
.collect();
let surf_refs: Vec<&Surface> = surfaces.iter().collect();
let (dx, dy) = spacing(self.opts.area_m2, self.ni, self.nj);
let stacked = layer_grid_stack(&surf_refs, dx, dy, &specs, self.spec.collapse_below_m)?;
let mut warnings = self.warnings.clone();
if stacked.truncated_cells > 0 {
warnings.push(BuildWarning::LayersTruncated {
cells: stacked.truncated_cells,
});
}
if stacked.collapsed_cells > 0 {
warnings.push(BuildWarning::CellsCollapsed {
cells: stacked.collapsed_cells,
});
}
if stacked.nk_capped {
warnings.push(BuildWarning::LayerCountCapped { nk: stacked.nk });
}
let mut grid = stacked.grid;
populate(
&mut grid,
self.opts.priors,
self.logs.as_deref(),
self.trend.as_ref(),
)?;
let mut property_reports = Vec::with_capacity(self.properties.len());
for pipe in &self.properties {
property_reports.push(pipe.apply_with_georef(&mut grid, self.spec.georef)?);
}
let zone_names: Vec<String> = stack.zone_layers.iter().map(|z| z.name.clone()).collect();
let zone_colors: Vec<Option<String>> =
stack.zone_layers.iter().map(|z| z.color.clone()).collect();
let zone_kranges: Vec<(String, core::ops::Range<usize>)> = zone_names
.iter()
.cloned()
.zip(
stacked
.zones
.iter()
.map(crate::gridder::StackedZone::k_range),
)
.collect();
self.apply_zone_population(&mut grid, &zone_kranges, &mut property_reports)?;
let horizon_names: Vec<String> = stack.horizons.iter().map(|h| h.name.clone()).collect();
let zone_contacts: Vec<Vec<Contact>> = stack
.zone_layers
.iter()
.map(|z| z.contacts.clone())
.collect();
let all_contacts: Vec<Contact> = zone_contacts.iter().flatten().copied().collect();
let ring = stack_boundary_ring(self.spec.boundary.as_ref(), self.spec.georef, nx, ny);
let framework = stack_framework(&stack.horizons, &surfaces, all_contacts, nx, ny, ring);
let per_zone_nk: Vec<usize> = stacked.zones.iter().map(|z| z.nk).collect();
let zones = ZoneTable::from_stack(&horizon_names, &zone_names, &zone_colors, &per_zone_nk);
let zone_prov: Vec<ZoneProvenance> = stacked
.zones
.iter()
.enumerate()
.map(|(z, sz)| ZoneProvenance {
name: zone_names[z].clone(),
top_horizon: horizon_names[z].clone(),
base_horizon: horizon_names[z + 1].clone(),
conformity: sz.conformity,
nk: sz.nk,
k_start: sz.k_start,
truncated_cells: sz.truncated_cells,
})
.collect();
let stack_prov = StackProvenance {
horizons: horizon_names,
zones: zone_prov,
interface_repairs,
};
let population = self.population_mode();
self.maybe_spill(
StaticModel::new(
framework,
grid,
zones,
Provenance {
inputs_ref: self.inputs_ref_string(),
solve_opts: self.opts.solve_opts,
conformity: self.opts.conformity,
nk: stacked.nk,
population,
realization: None,
warnings,
property_reports,
stack: Some(stack_prov),
well_ties,
sugar_cube: self.spec.sugar_cube,
},
self.spec.sw_gas,
)
.with_georef_opt(self.spec.georef)
.with_zone_contacts(Some(zone_contacts)),
)
}
fn apply_well_ties(
&self,
stack: &HorizonStack,
surfaces: &mut Vec<Surface>,
nx: usize,
ny: usize,
) -> Result<Vec<WellTieRecord>, StaticError> {
if self.spec.well_ties.is_empty() {
return Ok(Vec::new());
}
let (dx, dy) = spacing(self.opts.area_m2, self.ni, self.nj);
let mut tied_stack = stack.clone();
let (records, substituted) = substitute_tie_datums(
&mut tied_stack,
&self.spec.well_ties,
surfaces,
self.spec.ties,
(dx, dy),
nx,
ny,
)?;
if substituted {
*surfaces = resolve_stack_surfaces(&tied_stack, nx, ny, self.spec.extrapolation)?;
}
Ok(records)
}
}
pub(crate) fn substitute_tie_datums(
stack: &mut HorizonStack,
ties: &[WellTie],
surfaces: &[Surface],
settings: TieSettings,
(dx, dy): (f64, f64),
nx: usize,
ny: usize,
) -> Result<(Vec<WellTieRecord>, bool), StaticError> {
if let TieMethod::Radius { radius_m } = settings.method {
if !(radius_m.is_finite() && radius_m > 0.0) {
return Err(StaticError::InvalidInput(format!(
"tie radius must be finite and positive, got {radius_m}"
)));
}
}
let index_of = |name: &str, horizons: &[StackHorizon]| -> Option<usize> {
horizons.iter().position(|h| h.name == name)
};
let mut substituted = false;
let mut records = Vec::with_capacity(ties.len());
for tie in ties {
if tie.ip >= nx || tie.jp >= ny {
return Err(StaticError::Grid(format!(
"well tie '{}' node ({},{}) is off the {nx}x{ny} control lattice",
tie.id, tie.ip, tie.jp
)));
}
let mut residuals = Vec::with_capacity(tie.tops.len());
for (horizon, measured) in &tie.tops {
let h = index_of(horizon, &stack.horizons).ok_or_else(|| {
StaticError::InvalidInput(format!(
"well tie '{}' names unknown horizon '{horizon}'",
tie.id
))
})?;
let model_depth_m = surfaces[h].z(tie.ip, tie.jp);
let residual_m = *measured - model_depth_m;
residuals.push(HorizonTieResidual {
horizon: horizon.clone(),
measured_depth_m: *measured,
model_depth_m,
residual_m,
});
if let HorizonSource::Mapped(gd) = &mut stack.horizons[h].source {
let tie_idx = tie.jp * gd.ncol + tie.ip;
match settings.method {
TieMethod::Replace => {
gd.depth_m[tie_idx] = *measured;
gd.is_control[tie_idx] = true;
substituted = true;
}
TieMethod::Radius { radius_m } => {
for jp in 0..ny {
for ip in 0..nx {
let idx = jp * gd.ncol + ip;
if gd.depth_m[idx].is_nan() {
continue; }
let dist = ((ip as f64 - tie.ip as f64) * dx)
.hypot((jp as f64 - tie.jp as f64) * dy);
if dist <= radius_m {
let w = 1.0 - dist / radius_m;
gd.depth_m[idx] += residual_m * w;
}
}
}
gd.depth_m[tie_idx] = *measured;
gd.is_control[tie_idx] = true;
substituted = true;
}
}
}
}
records.push(WellTieRecord {
id: tie.id.clone(),
x: tie.x,
y: tie.y,
ip: tie.ip,
jp: tie.jp,
residuals,
});
}
Ok((records, substituted))
}
pub(crate) fn spacing(area_m2: f64, ni: usize, nj: usize) -> (f64, f64) {
let side = area_m2.sqrt();
(side / ni as f64, side / nj as f64)
}
pub(crate) fn condition_scatter(
stack: &mut HorizonStack,
frame: &StackFrame,
) -> Result<(), StaticError> {
let prof = std::env::var_os("SRS_PROFILE").is_some();
stack
.horizons
.par_iter_mut()
.try_for_each(|h| -> Result<(), StaticError> {
let HorizonSource::Scatter(points) = &h.source else {
return Ok(());
};
if points.is_empty() {
return Err(StaticError::InvalidInput(format!(
"scatter horizon '{}' has no points to grid",
h.name
)));
}
let t0 = std::time::Instant::now();
let gd = grid_scatter(points, frame);
if prof {
let ctrl = gd.is_control.iter().filter(|&&b| b).count();
eprintln!(
"[SRS_PROFILE] condition_scatter horizon='{}' points={} support_nodes={ctrl} ms={:.1}",
h.name,
points.len(),
t0.elapsed().as_secs_f64() * 1e3,
);
}
if gd.depth_m.iter().all(|z| z.is_nan()) {
return Err(StaticError::Grid(format!(
"scatter horizon '{}' landed no point on the {}x{} frame lattice \
(points outside the georeferenced extent)",
h.name,
frame.ni + 1,
frame.nj + 1
)));
}
h.source = HorizonSource::Mapped(gd);
Ok(())
})
}
pub(crate) struct ScatterConditioner {
op: MinCurvatureOperator,
support: Vec<bool>,
nx: usize,
ny: usize,
}
impl ScatterConditioner {
pub(crate) fn factor(points: &[WorldPoint], frame: &StackFrame) -> Option<(Self, Vec<f64>)> {
let (nx, ny) = (frame.ni + 1, frame.nj + 1);
let g = &frame.georef;
let mut sample_xy: Vec<[f64; 2]> = Vec::with_capacity(points.len());
let mut depths: Vec<f64> = Vec::with_capacity(points.len());
let mut support = vec![false; nx * ny];
for p in points {
if !(p.x.is_finite() && p.y.is_finite() && p.depth_m.is_finite()) {
continue;
}
let fi = (p.x - g.origin_x) / g.spacing_x + 0.5;
let fj = (p.y - g.origin_y) / g.spacing_y + 0.5;
if fi < -0.5 || fj < -0.5 || fi > nx as f64 - 0.5 || fj > ny as f64 - 0.5 {
continue; }
sample_xy.push([fi, fj]);
depths.push(p.depth_m);
let i0 = (fi.floor() as isize).clamp(0, nx as isize - 1) as usize;
let j0 = (fj.floor() as isize).clamp(0, ny as isize - 1) as usize;
for jp in j0..=(j0 + 1).min(ny - 1) {
for ip in i0..=(i0 + 1).min(nx - 1) {
support[jp * nx + ip] = true;
}
}
}
if sample_xy.is_empty() {
return None;
}
let lattice = Lattice::regular(0.0, 0.0, 1.0, 1.0, nx, ny);
let op = MinCurvatureOperator::factor(&lattice, &sample_xy, Conditioning::Bilinear).ok()?;
Some((
ScatterConditioner {
op,
support,
nx,
ny,
},
depths,
))
}
pub(crate) fn resolve(&self, depths: &[f64]) -> Option<Vec<f64>> {
let field = self.op.solve(depths).ok()?;
let mut depth_m = vec![f64::NAN; self.nx * self.ny];
for jp in 0..self.ny {
for ip in 0..self.nx {
if self.support[jp * self.nx + ip] {
depth_m[jp * self.nx + ip] = field[[ip, jp]];
}
}
}
Some(depth_m)
}
}
fn grid_scatter(points: &[WorldPoint], frame: &StackFrame) -> GriddedDepth {
let (nx, ny) = (frame.ni + 1, frame.nj + 1);
let (depth_m, support) = match ScatterConditioner::factor(points, frame) {
Some((cond, depths)) => match cond.resolve(&depths) {
Some(field) => (field, cond.support),
None => grid_scatter_sor_fallback(points, frame, cond.support),
},
None => (vec![f64::NAN; nx * ny], vec![false; nx * ny]),
};
GriddedDepth {
ncol: nx,
nrow: ny,
depth_m,
is_control: support,
}
}
fn grid_scatter_sor_fallback(
points: &[WorldPoint],
frame: &StackFrame,
support: Vec<bool>,
) -> (Vec<f64>, Vec<bool>) {
let (nx, ny) = (frame.ni + 1, frame.nj + 1);
let g = &frame.georef;
let mut coords: Vec<[f64; 3]> = Vec::with_capacity(points.len());
for p in points {
if !(p.x.is_finite() && p.y.is_finite() && p.depth_m.is_finite()) {
continue;
}
let fi = (p.x - g.origin_x) / g.spacing_x + 0.5;
let fj = (p.y - g.origin_y) / g.spacing_y + 0.5;
if fi < -0.5 || fj < -0.5 || fi > nx as f64 - 0.5 || fj > ny as f64 - 0.5 {
continue;
}
coords.push([fi, fj, p.depth_m]);
}
let mut depth_m = vec![f64::NAN; nx * ny];
if !coords.is_empty() {
let lattice = Lattice::regular(0.0, 0.0, 1.0, 1.0, nx, ny);
if let Ok(field) =
grid_min_curvature_conditioned(&coords, &lattice, None, Conditioning::Bilinear)
{
for jp in 0..ny {
for ip in 0..nx {
if support[jp * nx + ip] {
depth_m[jp * nx + ip] = field[[ip, jp]];
}
}
}
}
}
(depth_m, support)
}
pub(crate) fn validate_stack(stack: &HorizonStack) -> Result<(usize, usize), StaticError> {
let n = stack.horizons.len();
if n < 2 {
return Err(StaticError::InvalidInput(format!(
"a horizon stack needs at least 2 horizons, got {n}"
)));
}
if stack.zone_layers.len() != n - 1 {
return Err(StaticError::InvalidInput(format!(
"a {n}-horizon stack needs {} zones, got {}",
n - 1,
stack.zone_layers.len()
)));
}
let (nx, ny) = match &stack.horizons[0].source {
HorizonSource::Mapped(gd) => (gd.ncol, gd.nrow),
HorizonSource::TopsOnly(_) => {
return Err(StaticError::InvalidInput(
"the first (top) horizon must be a mapped surface".into(),
))
}
HorizonSource::Scatter(_) => {
return Err(StaticError::InvalidInput(
"raw-scatter horizons must be built through `from_scatter_stack` \
(they are conditioned to the model lattice before resolution)"
.into(),
))
}
};
if nx < 2 || ny < 2 {
return Err(StaticError::Grid(format!(
"top surface needs at least 2x2 nodes, got {nx}x{ny}"
)));
}
let mut seen_mapped = false;
for h in &stack.horizons {
match &h.source {
HorizonSource::Scatter(_) => {
return Err(StaticError::InvalidInput(format!(
"raw-scatter horizon '{}' must be conditioned via `from_scatter_stack`",
h.name
)))
}
HorizonSource::Mapped(gd) => {
if gd.ncol != nx || gd.nrow != ny {
return Err(StaticError::InvalidInput(format!(
"mapped horizon '{}' lattice {}x{} does not match the top {}x{}",
h.name, gd.ncol, gd.nrow, nx, ny
)));
}
if gd.depth_m.iter().all(|z| z.is_nan()) {
return Err(StaticError::Grid(format!(
"mapped horizon '{}' is fully undefined (all NaN)",
h.name
)));
}
seen_mapped = true;
}
HorizonSource::TopsOnly(picks) => {
if !seen_mapped {
return Err(StaticError::InvalidInput(format!(
"tops-only horizon '{}' has no mapped horizon above to drape from",
h.name
)));
}
if picks.is_empty() {
return Err(StaticError::InvalidInput(format!(
"tops-only horizon '{}' needs at least one pick",
h.name
)));
}
for p in picks {
if p.ip >= nx || p.jp >= ny {
return Err(StaticError::Grid(format!(
"tops-only horizon '{}' pick ({},{}) is off the {nx}x{ny} lattice",
h.name, p.ip, p.jp
)));
}
}
}
}
}
Ok((nx - 1, ny - 1))
}
pub(crate) fn pick_thickness_field(
picks: &[Pick],
above: &Surface,
nx: usize,
ny: usize,
policy: ExtrapolationPolicy,
) -> Result<Vec<f64>, StaticError> {
if picks.len() == 1 {
let p = picks[0];
let t = (p.depth_m - above.z(p.ip, p.jp)).max(0.0);
return Ok(vec![t; nx * ny]);
}
let controls: Vec<Control> = picks
.iter()
.map(|p| Control {
ip: p.ip,
jp: p.jp,
z: (p.depth_m - above.z(p.ip, p.jp)).max(0.0),
})
.collect();
let field_surf = warm_surface(nx, ny, &controls)?
.surface()
.taper_beyond_data(&controls, policy);
let mut field = vec![0.0; nx * ny];
for jp in 0..ny {
for ip in 0..nx {
field[jp * nx + ip] = field_surf.z(ip, jp).max(0.0);
}
}
Ok(field)
}
pub(crate) fn drape_tops_only(
picks: &[Pick],
above: &Surface,
nx: usize,
ny: usize,
policy: ExtrapolationPolicy,
) -> Result<Surface, StaticError> {
let field = pick_thickness_field(picks, above, nx, ny, policy)?;
above.offset_by_field(&field)
}
fn colocated_thickness_samples(
gd_a: &GriddedDepth,
gd_b: &GriddedDepth,
offset: Option<&[f64]>,
nx: usize,
ny: usize,
) -> Vec<Control> {
let mut samples = Vec::new();
for jp in 0..ny {
for ip in 0..nx {
let idx = jp * nx + ip;
let (za, zb) = (gd_a.depth_m[idx], gd_b.depth_m[idx]);
if !za.is_nan() && !zb.is_nan() {
let base = offset.map_or(0.0, |o| o[idx]);
samples.push(Control {
ip,
jp,
z: base + (zb - za).max(0.0),
});
}
}
}
samples
}
#[allow(clippy::needless_range_loop)] pub(crate) fn resolve_stack_surfaces(
stack: &HorizonStack,
nx: usize,
ny: usize,
policy: ExtrapolationPolicy,
) -> Result<Vec<Surface>, StaticError> {
let n = stack.horizons.len();
let mapped_idx: Vec<usize> = (0..n)
.filter(|&i| matches!(stack.horizons[i].source, HorizonSource::Mapped(_)))
.collect();
let top_i = mapped_idx[0];
let HorizonSource::Mapped(gd_top) = &stack.horizons[top_i].source else {
unreachable!("mapped_idx only indexes Mapped horizons");
};
let top_controls = controls_from_surface(gd_top);
if top_controls.is_empty() {
return Err(StaticError::Grid(format!(
"mapped horizon '{}' is fully undefined (all NaN)",
stack.horizons[top_i].name
)));
}
let top = warm_surface(nx, ny, &top_controls)?
.surface()
.taper_beyond_data(&top_controls, policy);
let mut surfaces: Vec<Option<Surface>> = vec![None; n];
let mut cums: Vec<Option<Vec<f64>>> = vec![None; n];
surfaces[top_i] = Some(top.clone());
cums[top_i] = Some(vec![0.0; nx * ny]);
let windows: Vec<[usize; 2]> = mapped_idx.windows(2).map(|w| [w[0], w[1]]).collect();
let primary_iso: Vec<Option<Vec<f64>>> = windows
.par_iter()
.map(|&[_, b]| -> Result<Option<Vec<f64>>, StaticError> {
let HorizonSource::Mapped(gd_b) = &stack.horizons[b].source else {
unreachable!("mapped_idx only indexes Mapped horizons");
};
let samples = colocated_thickness_samples(gd_top, gd_b, None, nx, ny);
if samples.is_empty() {
return Ok(None); }
let solved = warm_surface(nx, ny, &samples)?
.surface()
.taper_beyond_data(&samples, policy);
Ok(Some(
(0..nx * ny)
.map(|idx| {
let (ip, jp) = (idx % nx, idx / nx);
solved.z(ip, jp).max(0.0)
})
.collect(),
))
})
.collect::<Result<Vec<_>, StaticError>>()?;
for (wi, w) in mapped_idx.windows(2).enumerate() {
let (a, b) = (w[0], w[1]);
let HorizonSource::Mapped(gd_b) = &stack.horizons[b].source else {
unreachable!("mapped_idx only indexes Mapped horizons");
};
let cum_a = cums[a].clone().expect("upper anchor cum built");
let iso: Vec<f64> = if let Some(pre) = primary_iso[wi].clone() {
pre
} else {
let mut samples: Vec<Control> = Vec::new();
for &a2 in mapped_idx.iter().filter(|&&i| i < b).rev() {
if a2 == top_i {
break; }
let HorizonSource::Mapped(gd_a2) = &stack.horizons[a2].source else {
unreachable!("mapped_idx only indexes Mapped horizons");
};
let cum_a2 = cums[a2].as_deref().expect("shallower cum built");
samples = colocated_thickness_samples(gd_a2, gd_b, Some(cum_a2), nx, ny);
if !samples.is_empty() {
break;
}
}
if samples.is_empty() {
let controls_b = controls_from_surface(gd_b);
if controls_b.is_empty() {
return Err(StaticError::Grid(format!(
"mapped horizon '{}' is fully undefined (all NaN)",
stack.horizons[b].name
)));
}
let solved = warm_surface(nx, ny, &controls_b)?
.surface()
.taper_beyond_data(&controls_b, policy);
let above = surfaces[a].as_ref().expect("upper anchor built");
(0..nx * ny)
.map(|idx| {
let (ip, jp) = (idx % nx, idx / nx);
(solved.z(ip, jp) - above.z(ip, jp)).max(0.0)
})
.collect()
} else {
let solved = warm_surface(nx, ny, &samples)?
.surface()
.taper_beyond_data(&samples, policy);
(0..nx * ny)
.map(|idx| {
let (ip, jp) = (idx % nx, idx / nx);
solved.z(ip, jp).max(0.0)
})
.collect()
}
};
let cum_b: Vec<f64> = cum_a.iter().zip(&iso).map(|(p, i)| i.max(*p)).collect();
surfaces[b] = Some(top.offset_by_field(&cum_b)?);
let above = surfaces[a].clone().expect("upper anchor built");
for k in (a + 1)..b {
let HorizonSource::TopsOnly(picks) = &stack.horizons[k].source else {
unreachable!("a horizon between two mapped anchors must be tops-only");
};
let split = pick_thickness_field(picks, &above, nx, ny, policy)?;
let cum_k: Vec<f64> = split
.iter()
.zip(cum_a.iter().zip(&cum_b))
.map(|(s, (ca, cb))| ca + s.min(cb - ca))
.collect();
surfaces[k] = Some(top.offset_by_field(&cum_k)?);
}
cums[b] = Some(cum_b);
}
let last_mapped = *mapped_idx
.last()
.expect("validate_stack requires a mapped top");
for k in (last_mapped + 1)..n {
let HorizonSource::TopsOnly(picks) = &stack.horizons[k].source else {
unreachable!("horizons below the last mapped anchor are tops-only");
};
let above = surfaces[last_mapped]
.clone()
.expect("last mapped anchor built");
surfaces[k] = Some(drape_tops_only(picks, &above, nx, ny, policy)?);
}
let surfaces: Vec<Surface> = surfaces
.into_iter()
.map(|s| s.expect("every horizon resolved"))
.collect();
Ok(surfaces)
}
#[allow(clippy::too_many_arguments)]
pub(crate) fn repair_interface(
upper: &Surface,
lower: &Surface,
min_thickness: Option<f64>,
clamp: bool,
upper_derived: bool,
lower_derived: bool,
interface: usize,
repairs: &mut Vec<InterfaceRepair>,
) -> Result<(Surface, Surface), StaticError> {
if upper_derived && !lower_derived {
match min_thickness {
Some(min_t) => {
let (rep, columns, worst_m) =
upper.repair_min_thickness_from_below(lower, min_t)?;
if columns > 0 {
repairs.push(InterfaceRepair {
interface,
columns,
worst_m,
});
}
Ok((rep, lower.clone()))
}
None => Ok((upper.guard_above(lower, clamp)?, lower.clone())),
}
} else {
match min_thickness {
Some(min_t) => {
let (rep, columns, worst_m) = lower.repair_min_thickness(upper, min_t)?;
if columns > 0 {
repairs.push(InterfaceRepair {
interface,
columns,
worst_m,
});
}
Ok((upper.clone(), rep))
}
None => Ok((upper.clone(), lower.guard_below(upper, clamp)?)),
}
}
}
pub(crate) fn stack_boundary_ring(
explicit: Option<&Vec<[f64; 2]>>,
georef: Option<Georef>,
nx: usize,
ny: usize,
) -> Vec<[f64; 2]> {
if let Some(ring) = explicit {
return ring.clone();
}
if let Some(g) = georef {
let (ni, nj) = ((nx - 1) as f64, (ny - 1) as f64);
let xmin = g.origin_x - g.spacing_x / 2.0;
let xmax = g.origin_x + (ni - 0.5) * g.spacing_x;
let ymin = g.origin_y - g.spacing_y / 2.0;
let ymax = g.origin_y + (nj - 0.5) * g.spacing_y;
return vec![
[xmin, ymin],
[xmax, ymin],
[xmax, ymax],
[xmin, ymax],
[xmin, ymin],
];
}
vec![[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]]
}
pub(crate) fn stack_framework(
horizons: &[StackHorizon],
surfaces: &[Surface],
contacts: Vec<Contact>,
nx: usize,
ny: usize,
ring: Vec<[f64; 2]>,
) -> Wireframe {
let names: Vec<String> = horizons.iter().map(|h| h.name.clone()).collect();
stack_framework_from_names(&names, surfaces, contacts, nx, ny, ring)
}
pub(crate) fn stack_framework_from_names(
names: &[String],
surfaces: &[Surface],
contacts: Vec<Contact>,
nx: usize,
ny: usize,
ring: Vec<[f64; 2]>,
) -> Wireframe {
let n = names.len();
let hs: Vec<Horizon> = names
.iter()
.enumerate()
.map(|(i, name)| {
let role = if i == 0 {
HorizonRole::Top
} else if i == n - 1 {
HorizonRole::Base
} else {
HorizonRole::Intermediate
};
let mut depth = vec![0.0; nx * ny];
for jp in 0..ny {
for ip in 0..nx {
depth[jp * nx + ip] = surfaces[i].z(ip, jp);
}
}
Horizon {
name: name.clone(),
role,
surface: GriddedDepth {
ncol: nx,
nrow: ny,
depth_m: depth,
is_control: vec![true; nx * ny],
},
}
})
.collect();
Wireframe {
boundary: crate::wireframe::Boundary {
ring,
hardness: Hardness::Interpolated,
},
horizons: std::sync::Arc::new(hs),
contacts,
}
}
pub(crate) fn warm_surface(
nx: usize,
ny: usize,
controls: &[Control],
) -> Result<KernelSurface, StaticError> {
if controls.is_empty() {
return Err(StaticError::InvalidInput(
"structural solve needs at least one control point".into(),
));
}
crate::gridder::solve_surface_converged(nx, ny, controls)
}
fn flat_corner_surface(
nx: usize,
ny: usize,
controls: &[Control],
) -> Option<Result<KernelSurface, StaticError>> {
if controls.len() != 4 {
return None;
}
let z = controls[0].z;
if !z.is_finite() || controls.iter().any(|c| c.z != z) {
return None;
}
let max_i = nx.checked_sub(1)?;
let max_j = ny.checked_sub(1)?;
let corners = [(0, 0), (max_i, 0), (0, max_j), (max_i, max_j)];
if corners
.iter()
.all(|&(ip, jp)| controls.iter().any(|c| c.ip == ip && c.jp == jp))
{
Some(Ok(KernelSurface::flat(nx, ny, z)))
} else {
None
}
}
pub(crate) struct WireframeStructure {
pub ni: usize,
pub nj: usize,
pub top_controls: Vec<Control>,
pub base_controls: Option<Vec<Control>>,
pub warnings: Vec<BuildWarning>,
}
pub(crate) fn controls_from_surface(s: &GriddedDepth) -> Vec<Control> {
let mut controls = Vec::new();
for r in 0..s.nrow {
for c in 0..s.ncol {
let z = s.depth_m[r * s.ncol + c];
if !z.is_nan() {
controls.push(Control { ip: c, jp: r, z });
}
}
}
controls
}
pub(crate) fn wireframe_structure(
wf: &Wireframe,
opts: &BuildOpts,
) -> Result<WireframeStructure, StaticError> {
let top = wf
.horizons
.iter()
.find(|h| h.role == HorizonRole::Top)
.ok_or_else(|| StaticError::Grid("wireframe has no Top horizon".into()))?;
let s = &top.surface;
if s.ncol < 2 || s.nrow < 2 {
return Err(StaticError::Grid(
"top surface needs at least 2x2 nodes".into(),
));
}
if opts.nk == 0 {
return Err(StaticError::Grid("nk must be non-zero".into()));
}
validate_positive("area_m2", opts.area_m2)?;
validate_positive("gross_height_m", opts.gross_height_m)?;
let top_controls = controls_from_surface(s);
if top_controls.is_empty() {
return Err(StaticError::Grid(
"top surface is fully undefined (all NaN)".into(),
));
}
let mut base_controls = None;
let mut warnings = Vec::new();
for h in wf.horizons.iter() {
match h.role {
HorizonRole::Top => {}
HorizonRole::Base => {
let bs = &h.surface;
let reason = if bs.ncol != s.ncol || bs.nrow != s.nrow {
Some(format!(
"Base lattice {}x{} does not match Top {}x{}",
bs.ncol, bs.nrow, s.ncol, s.nrow
))
} else if base_controls.is_some() {
Some(
"a Base horizon is already in use (only the first is honoured)".to_string(),
)
} else {
let bc = controls_from_surface(bs);
if bc.is_empty() {
Some("Base surface is fully undefined (all NaN)".to_string())
} else {
base_controls = Some(bc);
None
}
};
if let Some(reason) = reason {
warnings.push(BuildWarning::UnusedHorizon {
name: h.name.clone(),
role: h.role,
reason,
});
}
}
HorizonRole::Intermediate => warnings.push(BuildWarning::UnusedHorizon {
name: h.name.clone(),
role: h.role,
reason: "intermediate horizons are unused until multi-zone layering (P5)"
.to_string(),
}),
}
}
if wf.contacts.is_empty() {
return Err(StaticError::Grid("wireframe has no fluid contact".into()));
}
Ok(WireframeStructure {
ni: s.ncol - 1,
nj: s.nrow - 1,
top_controls,
base_controls,
warnings,
})
}
fn synth_framework(ncol: usize, nrow: usize, top_depth_m: f64, contact_depth_m: f64) -> Wireframe {
Wireframe {
boundary: crate::wireframe::Boundary {
ring: vec![[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]],
hardness: Hardness::Assumed,
},
horizons: std::sync::Arc::new(vec![Horizon {
name: "Top".to_string(),
role: HorizonRole::Top,
surface: crate::wireframe::GriddedDepth {
ncol,
nrow,
depth_m: vec![top_depth_m; ncol * nrow],
is_control: vec![false; ncol * nrow],
},
}]),
contacts: vec![Contact {
kind: ContactKind::Owc,
depth_m: contact_depth_m,
hardness: Hardness::Assumed,
}],
}
}
#[cfg(test)]
mod conformity_tests {
use super::*;
use crate::gridder::Conformity;
use crate::model::draw::RealizationDraw;
use crate::model::{
Gaussian, PropertyPipeline, SectionSpec, StaticModelTemplate, UpscaleMethod, WellLog,
};
use crate::wireframe::{
Boundary, Contact, ContactKind, GriddedDepth, Hardness, Horizon, HorizonRole, Wireframe,
};
use petektools::{Variogram, VariogramModel};
const DZ: f64 = 20.0;
fn wedge_wf(n: usize, contacts: Vec<Contact>) -> Wireframe {
let top = vec![5000.0; n * n];
let mut base = vec![0.0; n * n];
for r in 0..n {
for c in 0..n {
base[r * n + c] = 5000.0 + 20.0 + 180.0 * (c as f64 / (n as f64 - 1.0));
}
}
let surf = |depth_m: Vec<f64>| GriddedDepth {
ncol: n,
nrow: n,
depth_m,
is_control: vec![true; n * n],
};
Wireframe {
boundary: Boundary {
ring: vec![[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]],
hardness: Hardness::Hard,
},
horizons: std::sync::Arc::new(vec![
Horizon {
name: "TopRes".into(),
role: HorizonRole::Top,
surface: surf(top),
},
Horizon {
name: "BaseRes".into(),
role: HorizonRole::Base,
surface: surf(base),
},
]),
contacts,
}
}
fn owc(depth_m: f64) -> Vec<Contact> {
vec![Contact {
kind: ContactKind::Owc,
depth_m,
hardness: Hardness::Hard,
}]
}
fn bopts(conformity: Conformity) -> BuildOpts {
BuildOpts {
area_m2: 90_000.0, gross_height_m: 100.0,
nk: 8, conformity,
solve_opts: SolveOpts::default(),
priors: ConstantPriors {
porosity: 0.25,
net_to_gross: 0.8,
water_saturation: 0.3,
},
}
}
fn truncated_count(m: &StaticModel) -> usize {
m.provenance()
.warnings
.iter()
.find_map(|w| match w {
BuildWarning::LayersTruncated { cells } => Some(*cells),
_ => None,
})
.unwrap_or(0)
}
#[test]
fn follow_styles_conserve_grv_and_hcpv_vs_proportional() {
let ip = |c| {
StaticModelBuilder::from_wireframe(&wedge_wf(11, owc(9000.0)), bopts(c))
.unwrap()
.build()
.unwrap()
.in_place()
.unwrap()
};
let prop = ip(Conformity::Proportional);
let ftop = ip(Conformity::FollowTop { dz_m: DZ });
let fbase = ip(Conformity::FollowBase { dz_m: DZ });
for (name, r) in [("FollowTop", &ftop), ("FollowBase", &fbase)] {
assert!(
(r.grv_m3 - prop.grv_m3).abs() / prop.grv_m3 < 1e-9,
"{name} GRV {} != proportional {}",
r.grv_m3,
prop.grv_m3
);
assert!(
(r.hcpv_m3 - prop.hcpv_m3).abs() / prop.hcpv_m3 < 1e-9,
"{name} HCPV {} != proportional {}",
r.hcpv_m3,
prop.hcpv_m3
);
}
}
#[test]
fn follow_top_truncates_and_populates_active_cells_without_nan() {
let vgm = Variogram::new(VariogramModel::Spherical, 0.0, 1.0, 120.0).unwrap();
let well_lo = WellLog::new(15.0, 15.0, vec![(5010.0, 0.18), (5060.0, 0.20)]);
let well_hi = WellLog::new(285.0, 285.0, vec![(5010.0, 0.24), (5180.0, 0.26)]);
let pipe = PropertyPipeline::new(crate::volumetrics::PORO)
.upscale(vec![well_lo, well_hi], UpscaleMethod::Arithmetic)
.propagate(Gaussian::new(vgm, 7).allow_mean_fill());
let m = StaticModelBuilder::from_wireframe(
&wedge_wf(11, owc(9000.0)),
bopts(Conformity::FollowTop { dz_m: DZ }),
)
.unwrap()
.with_property(pipe)
.build()
.unwrap();
assert_eq!(m.provenance().nk, 10, "dz-derived nk = ceil(200/20)");
assert!(truncated_count(&m) > 0, "thin updip columns truncate");
let grid = m.grid();
let dims = grid.dims();
let (ni, nj) = (dims.ni, dims.nj);
let cube = |n: &str| m.property(n).unwrap().values.as_slice();
let (poro, ntg, sw) = (
cube(crate::volumetrics::PORO),
cube(crate::volumetrics::NTG),
cube(crate::volumetrics::SW),
);
let mut active = 0usize;
for c in dims.iter() {
if grid.cell(c).dz() <= 1e-9 {
continue; }
active += 1;
let idx = (c.k * nj + c.j) * ni + c.i;
assert!(poro[idx].is_finite(), "PORO NaN in active cell {c:?}");
assert!(ntg[idx].is_finite(), "NTG NaN in active cell {c:?}");
assert!(sw[idx].is_finite(), "SW NaN in active cell {c:?}");
}
assert!(active > 0);
}
#[test]
fn two_contact_split_on_truncated_columns() {
let contacts = vec![
Contact {
kind: ContactKind::Goc,
depth_m: 5050.0,
hardness: Hardness::Hard,
},
Contact {
kind: ContactKind::Owc,
depth_m: 5120.0,
hardness: Hardness::Hard,
},
];
let m = StaticModelBuilder::from_wireframe(
&wedge_wf(11, contacts),
bopts(Conformity::FollowTop { dz_m: DZ }),
)
.unwrap()
.build()
.unwrap();
assert!(truncated_count(&m) > 0);
let ip = m.in_place().unwrap();
let gas = ip.gas.expect("gas cap");
let oil = ip.oil.expect("oil leg");
assert!(gas.hcpv_m3.is_finite() && oil.hcpv_m3.is_finite());
assert!(gas.grv_m3 > 0.0 && oil.grv_m3 > 0.0);
assert!(
(ip.hcpv_m3 - (gas.hcpv_m3 + oil.hcpv_m3)).abs() < 1e-6,
"total HCPV = gas + oil"
);
assert!(ip.per_cell_hcpv.iter().all(|v| v.is_finite()));
}
#[test]
fn template_realize_is_deterministic_under_follow_styles() {
let draw = RealizationDraw::new(90_000.0, 100.0, 9000.0, 0.25, 0.8, 0.3, 5);
for style in [
Conformity::FollowTop { dz_m: DZ },
Conformity::FollowBase { dz_m: DZ },
] {
let wf = wedge_wf(11, owc(9000.0));
let mut t1 = StaticModelTemplate::new(&wf, bopts(style)).unwrap();
let mut t2 = StaticModelTemplate::new(&wf, bopts(style)).unwrap();
let a = t1.realize(&draw).unwrap();
let b = t2.realize(&draw).unwrap();
assert_eq!(a.provenance().nk, 10, "dz-derived nk");
assert_eq!(a.grid().bulk_volume(), b.grid().bulk_volume());
assert_eq!(truncated_count(&a), truncated_count(&b));
assert!(
truncated_count(&a) > 0,
"{style:?} truncates the thin columns"
);
let (ia, ib) = (a.in_place().unwrap(), b.in_place().unwrap());
assert_eq!(ia.grv_m3, ib.grv_m3);
assert_eq!(ia.hcpv_m3, ib.hcpv_m3);
}
}
#[test]
fn section_bundle_nan_marks_inactive_layers() {
let m = StaticModelBuilder::from_wireframe(
&wedge_wf(11, owc(9000.0)),
bopts(Conformity::FollowTop { dz_m: DZ }),
)
.unwrap()
.build()
.unwrap();
let nk = m.provenance().nk;
let spec = SectionSpec::Polyline(vec![[15.0, 15.0], [285.0, 15.0]]);
let b = m
.intersection_bundle(&spec, Some(crate::volumetrics::PORO))
.unwrap();
assert!(!b.columns.is_empty());
let mut saw_nan = false;
for col in &b.columns {
assert_eq!(col.layer_tops.len(), nk, "arrays stay nk-sized");
assert_eq!(col.layer_bases.len(), nk);
assert_eq!(col.values.len(), nk);
for k in 0..nk {
let t = col.layer_tops[k];
if t.is_nan() {
saw_nan = true;
assert!(col.layer_bases[k].is_nan(), "base NaN with top");
assert!(col.values[k].is_nan(), "value NaN with geometry");
} else {
assert!(col.layer_bases[k].is_finite());
assert!(col.layer_bases[k] >= t);
assert!(col.values[k].is_finite(), "active PORO finite");
}
}
}
assert!(saw_nan, "the truncated updip columns must emit NaN layers");
}
#[test]
fn section_bundle_inactive_depths_serialize_to_json_null_not_zero() {
let m = StaticModelBuilder::from_wireframe(
&wedge_wf(11, owc(9000.0)),
bopts(Conformity::FollowTop { dz_m: DZ }),
)
.unwrap()
.build()
.unwrap();
let spec = SectionSpec::Polyline(vec![[15.0, 15.0], [285.0, 15.0]]);
let b = m
.intersection_bundle(&spec, Some(crate::volumetrics::PORO))
.unwrap();
let v = serde_json::to_value(&b).unwrap();
let cols = v.get("columns").and_then(|c| c.as_array()).unwrap();
let mut saw_null = false;
for (ci, col) in b.columns.iter().enumerate() {
let jtops = cols[ci]
.get("layer_tops")
.and_then(|t| t.as_array())
.unwrap();
for (k, &t) in col.layer_tops.iter().enumerate() {
if t.is_nan() {
saw_null = true;
assert!(
jtops[k].is_null(),
"inactive layer top must serialize to null"
);
} else {
assert!(
jtops[k].is_number(),
"active layer top must be a JSON number"
);
}
}
}
assert!(
saw_null,
"truncated columns must serialize inactive depths as null"
);
if let Some(traces) = v.get("horizon_traces").and_then(|t| t.as_array()) {
for tr in traces {
if let Some(depths) = tr.get("depths").and_then(|d| d.as_array()) {
for d in depths {
assert!(
d.is_null() || d.is_number(),
"a horizon-trace depth is either a number or null (never 0-for-gap)"
);
}
}
}
}
}
}
#[cfg(test)]
mod repair_precedence_tests {
use super::*;
use crate::gridder::Surface;
fn surf(z: Vec<f64>) -> Surface {
warm_surface(
2,
2,
&[
Control {
ip: 0,
jp: 0,
z: z[0],
},
Control {
ip: 1,
jp: 0,
z: z[1],
},
Control {
ip: 0,
jp: 1,
z: z[2],
},
Control {
ip: 1,
jp: 1,
z: z[3],
},
],
)
.unwrap()
.surface()
.clone()
}
#[test]
fn derived_upper_yields_to_mapped_lower() {
let upper = surf(vec![12.0, 5.0, 5.0, 5.0]);
let lower = surf(vec![10.0; 4]);
let mut reps = Vec::new();
let (u, l) =
repair_interface(&upper, &lower, None, true, true, false, 0, &mut reps).unwrap();
assert_eq!(l.z(0, 0), 10.0, "mapped lower is authoritative — untouched");
assert_eq!(
u.z(0, 0),
10.0,
"derived upper clamped up to the mapped lower"
);
assert_eq!(u.z(1, 0), 5.0, "non-crossing node untouched");
}
#[test]
fn mapped_lower_yields_when_both_mapped() {
let upper = surf(vec![12.0, 5.0, 5.0, 5.0]);
let lower = surf(vec![10.0; 4]);
let mut reps = Vec::new();
let (u, l) =
repair_interface(&upper, &lower, None, true, false, false, 0, &mut reps).unwrap();
assert_eq!(u.z(0, 0), 12.0, "upper untouched when both mapped");
assert_eq!(
l.z(0, 0),
12.0,
"lower clamped down to the upper (lower yields)"
);
}
}
#[cfg(test)]
mod scatter_conditioner_tests {
use super::*;
fn world_frame() -> StackFrame {
StackFrame {
ni: 12,
nj: 12,
georef: Georef::new(431_000.0 + 0.5 * 50.0, 6_521_000.0 + 0.5 * 50.0, 50.0, 50.0)
.unwrap(),
}
}
fn scatter(frame: &StackFrame, depth_at: impl Fn(f64, f64) -> f64) -> Vec<WorldPoint> {
let g = &frame.georef;
let mut pts = Vec::new();
for j in 0..11 {
for i in 0..11 {
let fi = i as f64 + 0.37;
let fj = j as f64 + 0.61;
let x = g.origin_x + (fi - 0.5) * g.spacing_x;
let y = g.origin_y + (fj - 0.5) * g.spacing_y;
pts.push(WorldPoint {
x,
y,
depth_m: depth_at(fi, fj),
});
}
}
pts
}
#[test]
fn resolve_on_a_reused_factor_is_bit_identical_to_conditioning_from_scratch() {
let frame = world_frame();
let pts_a = scatter(&frame, |fi, fj| 2000.0 + 3.0 * fi - 1.5 * fj);
let pts_b = scatter(&frame, |fi, fj| {
2000.0 + 3.0 * fi - 1.5 * fj + 7.0 * (0.3 * fi).sin() - 4.0 * (0.2 * fj).cos()
});
let (cond, depths_a) = ScatterConditioner::factor(&pts_a, &frame).expect("factor A");
let (_, depths_b) = ScatterConditioner::factor(&pts_b, &frame).expect("factor B");
let reused_a = cond.resolve(&depths_a).expect("resolve A");
let reused_b = cond.resolve(&depths_b).expect("resolve B");
let fresh_a = grid_scatter(&pts_a, &frame);
let fresh_b = grid_scatter(&pts_b, &frame);
let bit_identical = |x: &[f64], y: &[f64]| {
x.len() == y.len() && x.iter().zip(y).all(|(a, b)| a.to_bits() == b.to_bits())
};
assert!(
bit_identical(&reused_a, &fresh_a.depth_m),
"resolve on the reused factor must be bit-identical to conditioning A from scratch"
);
assert!(
bit_identical(&reused_b, &fresh_b.depth_m),
"resolve on the reused factor with B's depths must equal conditioning B from scratch \
(the fixed geometry, varying-depth MC lever)"
);
assert!(
reused_a
.iter()
.zip(&reused_b)
.any(|(a, b)| a.is_finite() && b.is_finite() && (a - b).abs() > 1e-6),
"the two draws must produce different conditioned surfaces"
);
}
}