+
+
+
+
+
+
+
+
+# SPDX-FileCopyrightText: 2026 NORCE Research AS
+# SPDX-License-Identifier: GPL-3.0
+# pylint: disable=C0103,R0902
+
+"""Configuration and simulation-data models shared across plopm workflows.
+
+PlopmConfig stores command-line options and normalized runtime settings used to
+create summary plots, spatial maps, animations, and VTK output. SimData stores
+OPM file handles, grid dimensions, and cell data loaded for one simulation case.
+
+Both objects are mutable because CLI values are normalized and simulation data
+are populated progressively during processing.
+"""
+
+from dataclasses import dataclass, field
+
+import numpy as np
+from numpy.typing import NDArray
+from opm.io.ecl import EclFile as OpmFile
+from opm.io.ecl import EGrid as OpmGrid
+from opm.io.ecl import ERst as OpmRestart
+
+
+
+[docs]
+@dataclass(slots=True)
+class PlopmConfig:
+ """Options and runtime settings for a plopm operation.
+
+ Most list fields contain one value per variable, case, or subplot after
+ initialization. Values read from the CLI are normalized before plotting so
+ downstream functions can use consistent indexing.
+
+ Attributes
+ ----------
+ gif, csv, png, vtk
+ Whether GIF, CSV, PNG, or VTK output is active for the current run.
+ equal_aspect
+ Whether spatial maps use the same scale along both coordinate axes.
+ remove_duplicate_labels
+ Whether repeated axis labels are hidden in subplot layouts.
+ list_variables
+ Whether available INIT, UNRST, and summary variables are printed.
+ gif_loop
+ Whether generated GIF animations repeat after the final frame.
+ step_plot
+ Whether one-dimensional series are drawn as step plots.
+ global_range
+ Whether map limits and features are evaluated globally instead of only
+ within the selected slice.
+ rst_range
+ Whether PNG color limits are evaluated over the restart range.
+ sensor
+ Whether one-dimensional values are extracted at a grid-cell sensor.
+ layer
+ Whether one-dimensional values are extracted along a grid axis or layer.
+ csv_column_summary
+ Whether a one-dimensional series is read from CSV columns.
+ discrete
+ Whether the current spatial quantity uses discrete color categories.
+ fontsize
+ Base font size used in generated figures.
+ mask_threshold
+ Threshold applied to the selected mask variable.
+ gif_interval
+ Delay between GIF frames.
+ stress_coefficient
+ Vertical stress coefficient used for caprock-integrity quantities.
+ xscale, yscale
+ Factors converting grid coordinates to the requested spatial units.
+ ensemble
+ Ensemble mode controlling uncertainty bands and bounding members.
+ ncolors
+ Number of case-dependent styles used for summary plots.
+ color_log_ticks
+ Tick values used on logarithmic colorbars.
+ case_labels
+ User-provided case names used in legends and ensemble labels.
+ cases
+ Simulation-case paths grouped as requested by the CLI.
+ dual_grid
+ Per-variable flags enabling dual-porosity grid handling.
+ subplot_grid
+ Requested subplot rows and columns.
+ variables
+ Variables or variable expressions requested for processing.
+ filters
+ Property-filter expressions applied per case or variable.
+ title
+ Per-plot titles.
+ clim
+ Lower and upper color limits for spatial maps.
+ figsize
+ Figure width and height for each generated plot.
+ min_threshold, max_threshold
+ Limits outside which quantity values are hidden.
+ grid_edges
+ Per-map settings controlling cell-edge drawing.
+ colorbar_tick_count
+ Requested number of colorbar ticks.
+ legend_labels
+ Labels for cases, variables, or ensemble bounds.
+ hide_map_elements
+ Map components to omit, such as axes, labels, or colorbars.
+ time_units
+ Requested time unit for each one-dimensional plot.
+ scale_factor
+ Multipliers applied to plotted or exported quantity values.
+ axis_grid
+ Per-plot settings controlling the Matplotlib axis grid.
+ dpi
+ Output resolution for each generated figure.
+ colorbar_ticks
+ Explicit colorbar tick values.
+ legend_location
+ Per-plot legend placement.
+ vtk_format
+ VTK data type selected for each exported variable.
+ vtk_names
+ Variable names written to VTK cell-data arrays.
+ color_log
+ Flags selecting logarithmic color normalization.
+ rotation
+ Rotation angles applied to grid coordinates, in degrees.
+ filename
+ Output filenames normalized per requested plot.
+ translation
+ Coordinate translations applied after rotation.
+ restart
+ Selected OPM restart report steps.
+ aggregation
+ Aggregation method applied through a slice or selected cells.
+ distance
+ Distance method and target, such as a sensor or model boundaries.
+ histogram
+ Histogram settings, including the requested bins.
+ xlabel, ylabel
+ Per-plot axis labels.
+ xformat, yformat
+ Format strings used for axis tick labels.
+ xtick_count, ytick_count
+ Requested numbers of major ticks.
+ xlog, ylog
+ Flags selecting logarithmic axes.
+ xlim, ylim
+ Per-plot axis limits.
+ vsum
+ Summary-variable expressions prepared for plotting.
+ summary
+ Loaded or derived summary-series values.
+ time
+ Time coordinates associated with summary values.
+ wells, faults
+ Parsed feature locations used in spatial maps.
+ slice
+ Parsed i, j, and k selections used by all workflows.
+ csv_columns
+ CSV column settings retained in parsed per-plot form.
+ mass_vars
+ Supported component-mass quantities.
+ summary_mass
+ Summary vectors converted from standard volume to mass.
+ mass_fracs
+ Supported component mass-fraction quantities.
+ caprock_vars
+ Supported caprock-integrity quantities.
+ linewidth_values
+ Default line widths before per-variable normalization.
+ units
+ Display units associated with requested quantities.
+ cb_formats
+ Normalized numeric formats used for colorbar labels.
+ colormaps
+ Colormaps assigned to spatial variables.
+ disc_colormaps
+ Available colormaps suitable for discrete values.
+ linestyle, linewidth, colors
+ Normalized styles used by summary plots.
+ colors_default, linestyle_default
+ Default style sequences used when none are supplied.
+ colorbar_position
+ Relative position and size of an explicitly placed colorbar axis.
+ difference_input
+ Second case, folder, or file used to calculate differences.
+ colors_raw
+ Color specification received from the CLI before normalization.
+ output_dir
+ Directory in which generated files are written.
+ case
+ Primary case path used for file detection and classification.
+ fill_between_style
+ Colors and opacity values for ensemble uncertainty bands.
+ colorbar_format
+ Colorbar format specification received from the CLI.
+ fc
+ Figure or axes face color.
+ inactive_color
+ Color assigned to inactive grid cells.
+ mask_variable
+ Variable used to mask spatial-map values.
+ suptitle
+ Figure-level title shared by all subplots.
+ colorbar_label
+ User-provided colorbar label.
+ slice_mode
+ Mode used when retaining wells or faults in an aggregated slice.
+ xunits, yunits
+ Requested spatial unit codes for both coordinate axes.
+ xunit, yunit
+ Formatted spatial unit labels shown on the axes.
+ slices
+ Normalized half-open ranges used for spatial slice aggregation.
+ csv_cols
+ Normalized CSV column indices used for gridded CSV data.
+ """
+
+ # Output modes and processing switches
+ gif: bool = False
+ csv: bool = False
+ png: bool = False
+ vtk: bool = False
+ equal_aspect: bool = False
+ remove_duplicate_labels: bool = False
+ list_variables: bool = False
+ gif_loop: bool = False
+ step_plot: bool = False
+ global_range: bool = False
+ rst_range: bool = False # Evaluate color limits across restart steps
+ sensor: bool = False
+ layer: bool = False
+ csv_column_summary: bool = False
+ discrete: bool = True
+
+ # Scalar plot and animation settings
+ fontsize: float = 0.0
+ mask_threshold: float = 0.0
+ gif_interval: float = 0.0
+ stress_coefficient: float = 0.0
+ xscale: float = 1.0
+ yscale: float = 1.0
+ ensemble: int = 0 # 0: off; 1: band; 2: bounds; 3: both
+ ncolors: int = 1
+
+ # Input cases, variables, and normalized selections
+ color_log_ticks: list = field(default_factory=list)
+ case_labels: list = field(default_factory=list) # Before path expansion
+ cases: list = field(default_factory=list) # Nested groups of case stems
+ dual_grid: list = field(default_factory=list)
+ subplot_grid: list = field(default_factory=list)
+ variables: list = field(default_factory=list)
+ filters: list = field(default_factory=list)
+
+ # Figure, subplot, and axis settings
+ title: list = field(default_factory=list)
+ clim: list = field(default_factory=list)
+ figsize: list = field(default_factory=list)
+ min_threshold: list = field(default_factory=list)
+ max_threshold: list = field(default_factory=list)
+
+ # Color, line, and map styling
+ grid_edges: list = field(default_factory=list)
+ colorbar_tick_count: list = field(default_factory=list)
+ legend_labels: list = field(default_factory=list)
+ hide_map_elements: list = field(default_factory=list)
+ time_units: list = field(default_factory=list)
+ scale_factor: list = field(default_factory=list)
+ axis_grid: list = field(default_factory=list)
+ dpi: list = field(default_factory=list)
+ colorbar_ticks: list = field(default_factory=list)
+ legend_location: list = field(default_factory=list)
+
+ # Generated output names and VTK settings
+ vtk_format: list = field(default_factory=list)
+ vtk_names: list = field(default_factory=list)
+ color_log: list = field(default_factory=list)
+ rotation: list = field(default_factory=list)
+ filename: list = field(default_factory=list)
+ translation: list = field(default_factory=list)
+ restart: list = field(default_factory=list)
+ aggregation: list = field(default_factory=list)
+ distance: list = field(default_factory=list)
+ histogram: list = field(default_factory=list)
+ xlabel: list = field(default_factory=list)
+ xformat: list = field(default_factory=list)
+ xtick_count: list = field(default_factory=list)
+ xlog: list = field(default_factory=list)
+ xlim: list = field(default_factory=list)
+ ylabel: list = field(default_factory=list)
+ yformat: list = field(default_factory=list)
+ ytick_count: list = field(default_factory=list)
+ ylog: list = field(default_factory=list)
+ ylim: list = field(default_factory=list)
+
+ # Summary data, features, and derived quantities
+ vsum: list = field(default_factory=list)
+ summary: list = field(default_factory=list)
+ time: list = field(default_factory=list)
+ wells: list = field(default_factory=list)
+ faults: list = field(default_factory=list)
+ slice: list = field(default_factory=list) # Parsed i, j, k selections
+ csv_columns: list = field(default_factory=list)
+ mass_vars: list = field(default_factory=list)
+ summary_mass: list = field(default_factory=list)
+ mass_fracs: list = field(default_factory=list)
+ caprock_vars: list = field(default_factory=list)
+ linewidth_values: list = field(default_factory=list)
+ units: list = field(default_factory=list)
+ cb_formats: list = field(default_factory=list)
+ colormaps: list = field(default_factory=list)
+ disc_colormaps: list = field(default_factory=list)
+ linestyle: list = field(default_factory=list)
+ linewidth: list = field(default_factory=list)
+ colors: list = field(default_factory=list)
+ colors_default: list = field(default_factory=list)
+ linestyle_default: list = field(default_factory=list)
+ colorbar_position: tuple[float, float, float, float] = (-1.0, -1.0, -1.0, -1.0)
+
+ # String options and values derived during initialization
+ difference_input: str = ""
+ colors_raw: str = "" # Before normalization into colors or colormaps
+ output_dir: str = ""
+ case: str = ""
+ fill_between_style: str = ""
+ colorbar_format: str = ""
+ fc: str = ""
+ inactive_color: str = ""
+ mask_variable: str = ""
+ suptitle: str = ""
+ colorbar_label: str = ""
+ slice_mode: str = "" # min keeps intersections; max keeps exact positions
+ xunits: str = ""
+ yunits: str = ""
+ xunit: str = ""
+ yunit: str = ""
+
+
+
+
+[docs]
+@dataclass(slots=True)
+class SimData:
+ """OPM readers, grid properties, and selected report steps for one case.
+
+ Arrays in global cell order use the full ``nx * ny * nz`` grid. Arrays in
+ active-cell order follow the indexing used by INIT and UNRST properties.
+
+ Attributes
+ ----------
+ init, unrst, grid
+ OPM readers for static properties, restart properties, and grid geometry.
+ porv
+ Pore volume in global cell order; inactive cells are non-positive.
+ dx, dy, dz
+ Cell dimensions in active-cell order.
+ active_pv
+ Pore volume in active-cell order.
+ active_idx
+ Mapping from global cell indices to active-cell indices.
+ steps, times
+ Selected restart report steps and their simulation times.
+ ncells, nsteps
+ Total grid-cell and available report-step counts.
+ nx, ny, nz
+ Grid dimensions along the i, j, and k axes.
+ """
+
+ # OPM file readers and grid geometry
+ init: OpmFile = None
+ unrst: OpmRestart = None
+ grid: OpmGrid = None
+
+ # Cell properties and global-to-active mapping
+ porv: NDArray = field(default_factory=lambda: np.array([]))
+ dx: NDArray = field(default_factory=lambda: np.array([]))
+ dy: NDArray = field(default_factory=lambda: np.array([]))
+ dz: NDArray = field(default_factory=lambda: np.array([]))
+ active_pv: NDArray = field(default_factory=lambda: np.array([]))
+ active_idx: NDArray = field(default_factory=lambda: np.array([]))
+
+ # Selected report steps and simulation times
+ steps: list = field(default_factory=list)
+ times: list = field(default_factory=list)
+
+ # Grid and report-step dimensions
+ ncells: int = 0
+ nsteps: int = 0
+ nx: int = 0
+ ny: int = 0
+ nz: int = 0
+
+
+
+
+
+
+# SPDX-FileCopyrightText: 2024-2026 NORCE Research AS
+# SPDX-License-Identifier: GPL-3.0
+# pylint: disable=C0302,R1702,W0123,W1401,R0912,R0914,R0915
+
+"""Command-line entry point and top-level workflow coordination for plopm.
+
+plopm supports three output workflows for OPM Flow simulation results:
+
+* One-dimensional plots and CSV files can be generated from summary vectors,
+ tabulated functions, grid-cell sensors, layers, histograms, and distances.
+* Two-dimensional PNG figures and GIF animations can be generated for selected
+ grid slices, with optional aggregation, masking, differences, wells, and faults.
+* VTK time series can be generated by combining OPM grid geometry with selected
+ INIT and UNRST properties.
+
+This module parses and validates command-line arguments, builds the runtime
+configuration, selects the appropriate workflow, and reports the generated
+files. Data reading, numerical processing, plotting, and file generation are
+implemented in the utility modules.
+"""
+
+import argparse
+import re
+import shlex
+import shutil
+import subprocess
+import sys
+
+from plopm.utils.initialization import (
+ build_config,
+ init_maps,
+ init_summary,
+ is_summary,
+)
+from plopm.utils.terminal import (
+ PlopmHelpFormatter,
+ cli_error_value,
+ plopm_error,
+ plopm_info,
+ plopm_name,
+ plopm_success,
+ plopm_tip,
+ warn_deprecated_options,
+)
+from plopm.utils.write_oned import make_plots
+from plopm.utils.write_twod import make_maps
+from plopm.utils.write_vtk import make_vtks
+
+
+
+[docs]
+def main(argv: list[str] | None = None) -> None:
+ """Run the plopm command-line workflow.
+
+ The function parses and validates CLI arguments, builds the shared
+ configuration, and dispatches VTK export, one-dimensional plotting, or
+ two-dimensional map generation. It reports the generated files after the
+ selected workflow completes.
+
+ Parameters
+ ----------
+ argv : list[str], optional
+ Arguments to parse instead of ``sys.argv[1:]``. This is primarily used
+ by tests and programmatic callers.
+
+ """
+ cmdargs = _load_parser(argv)
+ _check_cmdargs(cmdargs)
+ cfg = build_config(cmdargs)
+ if cfg.vtk:
+ plopm_info("processing, please wait...")
+ generated_files = make_vtks(
+ cmdargs.flow_path,
+ cfg.cases,
+ cfg.output_dir,
+ cfg.filename,
+ cfg.restart,
+ cfg.variables,
+ cfg.vtk_format,
+ cfg.vtk_names,
+ cfg.gif,
+ cfg.vtk,
+ cfg.filters,
+ cfg.scale_factor,
+ cfg.mass_vars,
+ cfg.mass_vars + cfg.mass_fracs,
+ cfg.caprock_vars,
+ cfg.stress_coefficient,
+ cfg.filters,
+ )
+ else:
+ if shutil.which("latex") is None:
+ plopm_tip(
+ "install LaTeX for improved fonts and text formatting; "
+ f"see the {plopm_name()} documentation for installation instructions."
+ )
+ if is_summary(cfg):
+ plopm_info("processing, please wait...")
+ init_summary(cfg)
+ generated_files = make_plots(cfg)
+ else:
+ plopm_info("processing, please wait...")
+ init_maps(cfg)
+ generated_files = make_maps(cfg)
+ plopm_success(cfg.output_dir, generated_files)
+
+
+
+def _load_parser(argv: list[str] | None = None) -> argparse.Namespace:
+ """Create the CLI parser and parse plopm arguments.
+
+ Parameters
+ ----------
+ argv : list[str], optional
+ Arguments to parse instead of ``sys.argv[1:]``.
+
+ Returns
+ -------
+ argparse.Namespace
+ Parsed command-line arguments.
+
+ """
+
+ parser = argparse.ArgumentParser(
+ formatter_class=PlopmHelpFormatter,
+ description=(
+ "plopm: Simplified and flexible Python tool for quick visualization "
+ "of OPM Flow geological models. See the online documentation for "
+ "examples and detailed option descriptions: "
+ "https://cssr-tools.github.io/plopm/introduction.html#option-reference"
+ ),
+ )
+
+ # ------------------------------------------------------------------
+ # Input and data selection
+ # ------------------------------------------------------------------
+
+ inputs = parser.add_argument_group("Input and data selection")
+
+ inputs.add_argument(
+ "-i",
+ "--input",
+ type=str.strip,
+ default="SPE11B",
+ help=(
+ "Base name or full path of the input. Separate multiple inputs "
+ 'with spaces, e.g. "SPE11B /home/user/SPE11B_TUNED"'
+ ),
+ )
+ inputs.add_argument(
+ "-v",
+ "--variable",
+ type=str.strip,
+ default="poro,permx,permz,porv,fipnum,satnum",
+ help=(
+ "Variable specification(s) to plot, including standard variables, "
+ "special variables, and expressions. Separate variables with commas"
+ ),
+ )
+ inputs.add_argument(
+ "-r",
+ "--restart",
+ type=str.strip,
+ default="-1",
+ help=(
+ "Restart step(s): a single step, comma-separated steps, or "
+ 'start:end[:step], e.g. "-1", "0,3,10", or "5:505:250"'
+ ),
+ )
+ inputs.add_argument(
+ "-cc",
+ "--csv-columns",
+ "-csv",
+ "--csv",
+ type=str.strip,
+ default="",
+ help=(
+ "CSV column indices starting at 1. Use t,value for time series or "
+ "x,y,value for spatial maps; separate inputs with semicolons"
+ ),
+ )
+ inputs.add_argument(
+ "-fp",
+ "--flow-path",
+ "-p",
+ "--path",
+ type=str.strip,
+ default="flow",
+ help="Path or command for the Flow executable used for VTK grid generation",
+ )
+
+ # ------------------------------------------------------------------
+ # Output options
+ # ------------------------------------------------------------------
+
+ output = parser.add_argument_group("Output options")
+
+ output.add_argument(
+ "-m",
+ "--format",
+ "--mode",
+ type=str.strip,
+ choices=["png", "gif", "csv", "vtk"],
+ default="png",
+ help="Output format",
+ )
+ output.add_argument(
+ "-o",
+ "--output-dir",
+ "--output",
+ type=str.strip,
+ default=".",
+ help="Base name or full path of the output directory",
+ )
+ output.add_argument(
+ "-fn",
+ "--filename",
+ "-save",
+ "--save",
+ type=str.strip,
+ default="",
+ help="Output file name",
+ )
+
+ # ------------------------------------------------------------------
+ # Spatial and temporal selection
+ # ------------------------------------------------------------------
+
+ selection = parser.add_argument_group("Spatial and temporal selection")
+
+ selection.add_argument(
+ "-s",
+ "--slice",
+ "--slide",
+ type=str.strip,
+ default=",1,",
+ help=(
+ "Spatial selection in i,j,k form, e.g. "
+ '"10,," for a plane, ",,5:10" for a range, '
+ '":,5,7" for a line, or "2,4,9" for a cell over time'
+ ),
+ )
+ selection.add_argument(
+ "-tu",
+ "--time-units",
+ "-tunits",
+ "--tunits",
+ type=str.strip,
+ choices=["s", "m", "h", "d", "w", "y", "dates", "empty", "tstep"],
+ default="d",
+ help="Summary-plot x-axis time units",
+ )
+ selection.add_argument(
+ "-dist",
+ "--distance",
+ "-distance",
+ type=str.strip,
+ choices=["min,sensor", "max,sensor", "min,border", "max,border", ""],
+ default="",
+ help="Compute the minimum or maximum distance to a sensor or lateral border",
+ )
+
+ # ------------------------------------------------------------------
+ # Filtering, masking, and thresholds
+ # ------------------------------------------------------------------
+
+ filtering = parser.add_argument_group("Filtering, masking, and thresholds")
+
+ filtering.add_argument(
+ "-flt",
+ "--filters",
+ "-filter",
+ type=str.strip,
+ default="",
+ help=(
+ "Cell-selection conditions. Join conditions for one input with '&' "
+ "and separate filters for different inputs with commas"
+ ),
+ )
+ filtering.add_argument(
+ "-vmin",
+ "--min-threshold",
+ "--vmin",
+ type=str.strip,
+ default="",
+ help="Minimum threshold used to remove variable values",
+ )
+ filtering.add_argument(
+ "-vmax",
+ "--max-threshold",
+ "--vmax",
+ type=str.strip,
+ default="",
+ help="Maximum threshold used to remove variable values",
+ )
+ filtering.add_argument(
+ "-mv",
+ "--mask-variable",
+ "-mask",
+ "--mask",
+ type=str.strip,
+ default="",
+ help="Static variable used as the background of a 2D map",
+ )
+ filtering.add_argument(
+ "-mt",
+ "--mask-threshold",
+ "-maskthr",
+ "--maskthr",
+ type=str.strip,
+ default="1e-3",
+ help="Threshold applied to the mask variable",
+ )
+
+ # ------------------------------------------------------------------
+ # Computation and data transformation
+ # ------------------------------------------------------------------
+
+ computation = parser.add_argument_group("Computation and data transformation")
+
+ computation.add_argument(
+ "-agg",
+ "--aggregation",
+ "-how",
+ "--how",
+ type=str.strip,
+ default="",
+ help=(
+ "Aggregation or selection method for 2D slices and projections: "
+ "min, max, sum, mean, pvmean, harmonic, arithmetic, first, or last"
+ ),
+ )
+ computation.add_argument(
+ "-sf",
+ "--scale-factor",
+ "-a",
+ "--adjust",
+ type=str.strip,
+ default="1",
+ help=(
+ "Multiplicative scaling factor applied to variable values, "
+ "e.g. 1e-9 to display mass in Mt"
+ ),
+ )
+ computation.add_argument(
+ "-di",
+ "--difference-input",
+ "-diff",
+ "--diff",
+ type=str.strip,
+ default="",
+ help="Base name or full path of the input model to subtract",
+ )
+ computation.add_argument(
+ "-sc",
+ "--stress-coefficient",
+ "-stress",
+ "--stress",
+ type=str.strip,
+ default="0.134",
+ help=(
+ "Stress coefficient used to compute pressure limits for "
+ "limipres, overpres, and objepres"
+ ),
+ )
+ computation.add_argument(
+ "-dg",
+ "--dual-grid",
+ "-dual",
+ "--dual",
+ type=str.strip,
+ default="0",
+ help="Enable dual-grid processing using 0 or 1",
+ )
+
+ # ------------------------------------------------------------------
+ # Plot types and statistical representation
+ # ------------------------------------------------------------------
+
+ plot_types = parser.add_argument_group("Plot types and statistical representation")
+
+ plot_types.add_argument(
+ "-hist",
+ "--histogram",
+ "-histogram",
+ type=str.strip,
+ default="",
+ help=(
+ "Histogram bins and optional distribution, e.g. "
+ '"20", "20,norm", or "20,lognorm"'
+ ),
+ )
+ plot_types.add_argument(
+ "-ens",
+ "--ensemble",
+ "-ensemble",
+ type=str.strip,
+ choices=["0", "1", "2", "3"],
+ default="0",
+ help=(
+ "Ensemble plotting mode: 0 disables it, 1 plots mean and error "
+ "bands, 2 plots minimum, mean, and maximum, and 3 plots both"
+ ),
+ )
+ plot_types.add_argument(
+ "-fb",
+ "--fill-between-style",
+ "-bandprop",
+ "--bandprop",
+ type=str.strip,
+ default="",
+ help="Fill colors and alpha values for ensemble error bands",
+ )
+ plot_types.add_argument(
+ "-sp",
+ "--step-plot",
+ "-step",
+ "--step",
+ type=str.strip,
+ choices=["0", "1"],
+ default="0",
+ help="Use ax.step instead of ax.plot",
+ )
+
+ # ------------------------------------------------------------------
+ # Figure and subplot layout
+ # ------------------------------------------------------------------
+
+ layout = parser.add_argument_group("Figure and subplot layout")
+
+ layout.add_argument(
+ "-fs",
+ "--figsize",
+ "-d",
+ "--dimensions",
+ type=str.strip,
+ default="7,5",
+ help='Figure width and height in inches, e.g. "8,16"',
+ )
+ layout.add_argument(
+ "-sg",
+ "--subplot-grid",
+ "-subfigs",
+ "--subfigs",
+ type=str.strip,
+ default="",
+ help='Number of subplot rows and columns, e.g. "2,2"',
+ )
+ layout.add_argument(
+ "-cbp",
+ "--colorbar-position",
+ "-cbsfax",
+ "--cbsfax",
+ type=str.strip,
+ default="0.2,0.01,0.6,0.02",
+ help=(
+ "Global colorbar position and size as left,bottom,width,height; "
+ "use 'empty' to remove it"
+ ),
+ )
+ layout.add_argument(
+ "-rdl",
+ "--remove-duplicate-labels",
+ "-delax",
+ "--delax",
+ type=str.strip,
+ choices=["0", "1"],
+ default="0",
+ help="Remove duplicated axis labels in subplot layouts",
+ )
+
+ # ------------------------------------------------------------------
+ # Titles, labels, and legends
+ # ------------------------------------------------------------------
+
+ text = parser.add_argument_group("Titles, labels, and legends")
+
+ text.add_argument(
+ "-t",
+ "--title",
+ type=str.strip,
+ default="0",
+ help="Figure title; separate titles for multiple plots with two spaces",
+ )
+ text.add_argument(
+ "-st",
+ "--suptitle",
+ "-suptitle",
+ type=str.strip,
+ default="",
+ help="Title for a group of subplots; use 0 to remove it",
+ )
+ text.add_argument(
+ "-xl",
+ "--xlabel",
+ "-xlabel",
+ type=str.strip,
+ default="",
+ help="X-axis label; separate labels for multiple plots with two spaces",
+ )
+ text.add_argument(
+ "-yl",
+ "--ylabel",
+ "-ylabel",
+ type=str.strip,
+ default="",
+ help="Y-axis label; separate labels for multiple plots with two spaces",
+ )
+ text.add_argument(
+ "-cbl",
+ "--colorbar-label",
+ "-clabel",
+ "--clabel",
+ type=str.strip,
+ default="",
+ help="Colorbar label; separate labels for multiple plots with two spaces",
+ )
+ text.add_argument(
+ "-llb",
+ "--legend-labels",
+ "-labels",
+ "--labels",
+ type=str.strip,
+ default="",
+ help="Summary-plot legend labels separated by two spaces",
+ )
+ text.add_argument(
+ "-ll",
+ "--legend-location",
+ "-loc",
+ "--loc",
+ type=str.strip,
+ default="best",
+ help="Legend location passed to matplotlib; use 'empty' to remove it",
+ )
+ text.add_argument(
+ "-hide",
+ "--hide-map-elements",
+ "-remove",
+ "--remove",
+ type=str.strip,
+ default="0,0,0,0",
+ help=(
+ "Hide the left axis, bottom axis, colorbar, and title using four "
+ "comma-separated values of 0 or 1"
+ ),
+ )
+
+ # ------------------------------------------------------------------
+ # Axes, coordinates, and formatting
+ # ------------------------------------------------------------------
+
+ axes = parser.add_argument_group("Axes, coordinates, and formatting")
+
+ axes.add_argument(
+ "-x",
+ "--xlim",
+ type=str.strip,
+ default="",
+ help='X-axis limits in display order, e.g. "[-100,200]"',
+ )
+ axes.add_argument(
+ "-y",
+ "--ylim",
+ type=str.strip,
+ default="",
+ help='Y-axis limits in display order, e.g. "[0,300]"',
+ )
+ axes.add_argument(
+ "-xu",
+ "--xunits",
+ "-xunits",
+ type=str.strip,
+ choices=["mm", "cm", "m", "km"],
+ default="m",
+ help="Spatial-map x-axis units",
+ )
+ axes.add_argument(
+ "-yu",
+ "--yunits",
+ "-yunits",
+ type=str.strip,
+ choices=["mm", "cm", "m", "km"],
+ default="m",
+ help="Spatial-map y-axis units",
+ )
+ axes.add_argument(
+ "-asp",
+ "--equal-aspect",
+ "-z",
+ "--scale",
+ type=str.strip,
+ choices=["0", "1"],
+ default="1",
+ help="Scale the axes equally in 2D maps",
+ )
+ axes.add_argument(
+ "-rot",
+ "--rotation",
+ "-rotate",
+ "--rotate",
+ type=str.strip,
+ default="0",
+ help="Grid rotation angle in degrees for 2D maps",
+ )
+ axes.add_argument(
+ "-tr",
+ "--translation",
+ "-translate",
+ "--translate",
+ type=str.strip,
+ default="[0,0]",
+ help='Grid translation in the x and y directions, e.g. "[100,-50]"',
+ )
+ axes.add_argument(
+ "-xlog",
+ "--xlog",
+ type=str.strip,
+ default="0",
+ help="Enable the logarithmic x-axis using 0 or 1",
+ )
+ axes.add_argument(
+ "-ylog",
+ "--ylog",
+ type=str.strip,
+ default="0",
+ help="Enable the logarithmic y-axis using 0 or 1",
+ )
+ axes.add_argument(
+ "-xf",
+ "--xformat",
+ "-xformat",
+ type=str.strip,
+ default="",
+ help='X-axis number format, e.g. ".2e"',
+ )
+ axes.add_argument(
+ "-yf",
+ "--yformat",
+ "-yformat",
+ type=str.strip,
+ default="",
+ help='Y-axis number format, e.g. ".1f"',
+ )
+ axes.add_argument(
+ "-xnt",
+ "--xtick-count",
+ "-xlnum",
+ "--xlnum",
+ type=str.strip,
+ default="5",
+ help="Number of x-axis ticks",
+ )
+ axes.add_argument(
+ "-ynt",
+ "--ytick-count",
+ "-ylnum",
+ "--ylnum",
+ type=str.strip,
+ default="5",
+ help="Number of y-axis ticks",
+ )
+
+ # ------------------------------------------------------------------
+ # Color scales and styling
+ # ------------------------------------------------------------------
+
+ styling = parser.add_argument_group("Color scales and styling")
+
+ styling.add_argument(
+ "-c",
+ "--colors",
+ type=str.strip,
+ default="",
+ help='Colormap or summary-plot colors, e.g. "jet" or "b,r"',
+ )
+ styling.add_argument(
+ "-cl",
+ "--clim",
+ "-b",
+ "--bounds",
+ type=str.strip,
+ default="",
+ help='Color-scale limits in display order, e.g. "[-0.1,11]"',
+ )
+ styling.add_argument(
+ "-clog",
+ "--color-log",
+ "-log",
+ "--log",
+ type=str.strip,
+ default="0",
+ help="Enable logarithmic color scaling using 0 or 1",
+ )
+ styling.add_argument(
+ "-clt",
+ "--color-log-ticks",
+ "-clogthks",
+ "--clogthks",
+ type=str.strip,
+ default="",
+ help='Tick values for logarithmic color scales, e.g. "[1,10,100]"',
+ )
+ styling.add_argument(
+ "-gr",
+ "--global-range",
+ "-global",
+ "--global",
+ type=str.strip,
+ choices=["0", "1"],
+ default="0",
+ help="Use the current slice range or whole-model range for color scaling",
+ )
+ styling.add_argument(
+ "-cbf",
+ "--colorbar-format",
+ "-cformat",
+ "--cformat",
+ type=str.strip,
+ default="",
+ help='Colorbar number format, e.g. ".2f"',
+ )
+ styling.add_argument(
+ "-cbn",
+ "--colorbar-tick-count",
+ "-cnum",
+ "--cnum",
+ type=str.strip,
+ default="",
+ help="Number of colorbar ticks",
+ )
+ styling.add_argument(
+ "-cbt",
+ "--colorbar-ticks",
+ "-cticks",
+ "--cticks",
+ type=str.strip,
+ default="",
+ help='Custom colorbar tick labels, e.g. "[A,B,C]"',
+ )
+ styling.add_argument(
+ "-lw",
+ "--linewidth",
+ "--lw",
+ type=str.strip,
+ default="",
+ help="Line widths separated by commas",
+ )
+ styling.add_argument(
+ "-ls",
+ "--linestyle",
+ "-e",
+ type=str.strip,
+ default="",
+ help='Line styles separated by commas, e.g. "solid,dotted"',
+ )
+ styling.add_argument(
+ "-ag",
+ "--axis-grid",
+ "-axgrid",
+ "--axgrid",
+ type=str.strip,
+ choices=["0", "1"],
+ default="1",
+ help="Display the summary-plot axis grid",
+ )
+ styling.add_argument(
+ "-fc",
+ "--facecolor",
+ "-facecolor",
+ type=str.strip,
+ default="w",
+ help="Color outside the spatial map",
+ )
+ styling.add_argument(
+ "-ic",
+ "--inactive-color",
+ "-ncolor",
+ "--ncolor",
+ type=str.strip,
+ default="w",
+ help="Color for inactive cells in 2D maps",
+ )
+ styling.add_argument(
+ "-ge",
+ "--grid-edges",
+ "-grid",
+ "--grid",
+ type=str.strip,
+ default="",
+ help="pcolormesh edge color and line width separated by a comma",
+ )
+ styling.add_argument(
+ "-fz",
+ "--fontsize",
+ "-f",
+ "--size",
+ type=str.strip,
+ default="12",
+ help="Font size",
+ )
+ styling.add_argument(
+ "-dpi",
+ "--dpi",
+ type=str.strip,
+ default="500",
+ help="Figure resolution in dots per inch",
+ )
+
+ # ------------------------------------------------------------------
+ # VTK output
+ # ------------------------------------------------------------------
+
+ vtk = parser.add_argument_group("VTK output")
+
+ vtk.add_argument(
+ "-vf",
+ "--vtk-format",
+ "-vtkformat",
+ "--vtkformat",
+ type=str.strip,
+ default="Float64",
+ help="VTK data type for each variable, separated by commas",
+ )
+ vtk.add_argument(
+ "-vn",
+ "--vtk-names",
+ "-vtknames",
+ "--vtknames",
+ type=str.strip,
+ default="",
+ help="Custom VTK variable names separated by commas",
+ )
+
+ # ------------------------------------------------------------------
+ # GIF output
+ # ------------------------------------------------------------------
+
+ gif = parser.add_argument_group("GIF output")
+
+ gif.add_argument(
+ "-gi",
+ "--gif-interval",
+ "-interval",
+ "--interval",
+ type=str.strip,
+ default="1000",
+ help="GIF frame interval in milliseconds",
+ )
+ gif.add_argument(
+ "-gl",
+ "--gif-loop",
+ "-loop",
+ "--loop",
+ type=str.strip,
+ default="0",
+ help="Loop GIF animations indefinitely using 0 or 1",
+ )
+
+ # ------------------------------------------------------------------
+ # Information and diagnostics
+ # ------------------------------------------------------------------
+
+ diagnostics = parser.add_argument_group("Information and diagnostics")
+
+ diagnostics.add_argument(
+ "-lv",
+ "--list-variables",
+ "-printv",
+ "--printv",
+ type=str.strip,
+ choices=["0", "1"],
+ default="0",
+ help="Print the available variables",
+ )
+
+ parsed_argv = sys.argv[1:] if argv is None else argv
+
+ warn_deprecated_options(parsed_argv)
+
+ return parser.parse_args(parsed_argv)
+
+
+def _check_cmdargs(cmdargs: argparse.Namespace) -> None:
+ """Validate command-line values and option combinations.
+
+ The function checks value syntax, accepted choices, mutually exclusive
+ operations, mode-specific options, and availability of OPM Flow for VTK
+ generation.
+
+ Parameters
+ ----------
+ cmdargs : argparse.Namespace
+ Parsed arguments returned by :func:`_load_parser`.
+
+ Raises
+ ------
+ SystemExit
+ If a value is invalid or incompatible options are requested.
+
+ """
+
+ mode = cmdargs.format
+ vtk_mode = mode == "vtk"
+ gif_mode = mode == "gif"
+ number = r"[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?"
+ positive_integer = r"[1-9]\d*"
+ non_negative_integer = r"\d+"
+
+ if not cmdargs.input:
+ plopm_error(f"the input {cli_error_value('-i')} cannot be empty.")
+ if not cmdargs.output_dir:
+ plopm_error(f"the output folder {cli_error_value('-o')} cannot be empty.")
+ if not cmdargs.variable:
+ plopm_error(f"the variable {cli_error_value('-v')} cannot be empty.")
+
+ positive_number_options = [
+ ("-fz", "fontsize"),
+ ("-dpi", "dpi"),
+ ("-xnt", "xtick_count"),
+ ("-ynt", "ytick_count"),
+ ("-mt", "mask_threshold"),
+ ("-gi", "gif_interval"),
+ ]
+ for option, name in positive_number_options:
+ raw_value = getattr(cmdargs, name)
+ value = _parse_number(option, raw_value)
+ if value <= 0:
+ plopm_error(
+ f"expected a positive number, not "
+ f"{cli_error_value(f'{option} {raw_value}')}."
+ )
+
+ number_options = [
+ ("-sc", "stress_coefficient"),
+ ("-rot", "rotation"),
+ ]
+ for option, name in number_options:
+ _parse_number(option, getattr(cmdargs, name))
+
+ __parse_number_list("-sf", cmdargs.scale_factor)
+
+ optional_number_options = [
+ ("-vmin", "min_threshold"),
+ ("-vmax", "max_threshold"),
+ ]
+ for option, name in optional_number_options:
+ value = getattr(cmdargs, name)
+ if value:
+ _parse_number(option, value)
+
+ if (
+ cmdargs.min_threshold
+ and cmdargs.max_threshold
+ and float(cmdargs.min_threshold) > float(cmdargs.max_threshold)
+ ):
+ plopm_error(
+ f"the minimum threshold "
+ f"{cli_error_value(f'-vmin {cmdargs.min_threshold}')} must not be "
+ f"greater than the maximum threshold "
+ f"{cli_error_value(f'-vmax {cmdargs.max_threshold}')}."
+ )
+
+ colorbar_tick_numbers = cmdargs.colorbar_tick_count
+ if colorbar_tick_numbers:
+ cnum_entries = colorbar_tick_numbers.split(",")
+ if any(not re.fullmatch(positive_integer, entry) for entry in cnum_entries):
+ plopm_error(
+ "expected positive integers separated by commas, "
+ f"not {cli_error_value(f'-cbn {colorbar_tick_numbers}')}."
+ )
+
+ boolean_options = [
+ ("-xlog", "xlog"),
+ ("-ylog", "ylog"),
+ ("-clog", "color_log"),
+ ("-dg", "dual_grid"),
+ ("-gl", "gif_loop"),
+ ]
+ for option, name in boolean_options:
+ raw_value = getattr(cmdargs, name)
+ values = raw_value.split(",")
+ if any(value not in ["0", "1"] for value in values):
+ plopm_error(
+ "expected values containing only 0 or 1, separated by commas, "
+ f"not {cli_error_value(f'{option} {raw_value}')}."
+ )
+
+ dimensions = __parse_number_list(
+ "-fs",
+ cmdargs.figsize,
+ 2,
+ )
+ if any(value <= 0 for value in dimensions):
+ plopm_error(
+ f"figure dimensions must be positive, not "
+ f"{cli_error_value(f'-fs {cmdargs.figsize}')}."
+ )
+
+ translation = cmdargs.translation
+ if not re.fullmatch(
+ rf"\[\s*{number}\s*,\s*{number}\s*\]",
+ translation,
+ ):
+ plopm_error(
+ f"expected two numbers enclosed by brackets, such as "
+ f"{cli_error_value('-tr [10,-5]')}, not "
+ f"{cli_error_value(f'-tr {translation}')}."
+ )
+
+ interval_pattern = re.compile(rf"\[\s*({number})\s*,\s*({number})\s*\]")
+ for option, name in [
+ ("-cl", "clim"),
+ ("-x", "xlim"),
+ ("-y", "ylim"),
+ ]:
+ value = getattr(cmdargs, name)
+ if not value:
+ continue
+ for interval_value in value.split():
+ if not interval_pattern.fullmatch(interval_value):
+ plopm_error(
+ f"expected two numeric bounds enclosed by brackets, such "
+ f"as {cli_error_value(f'{option} [0,10]')}, not "
+ f"{cli_error_value(f'{option} {interval_value}')}."
+ )
+
+ aggregation_methods = cmdargs.aggregation
+ if aggregation_methods:
+ valid_aggregation_methods = [
+ "min",
+ "max",
+ "sum",
+ "mean",
+ "pvmean",
+ "harmonic",
+ "arithmetic",
+ "first",
+ "last",
+ ]
+ method_entries = aggregation_methods.split(",")
+ if any(method not in valid_aggregation_methods for method in method_entries):
+ plopm_error(
+ f"expected aggregation methods from "
+ f"{', '.join(valid_aggregation_methods)}, not "
+ f"{cli_error_value(f'-agg {aggregation_methods}')}."
+ )
+
+ slice_value = cmdargs.slice
+ slices = slice_value.split()
+ slice_entry_pattern = re.compile(
+ rf"(?:{positive_integer}|" rf"{positive_integer}:{positive_integer}|:)?"
+ )
+ if not slices:
+ plopm_error(f"the slice selection {cli_error_value('-s')} cannot be empty.")
+
+ slice_entries: list[list[str]] = []
+ for selection in slices:
+ entries = selection.split(",")
+ if len(entries) != 3 or any(
+ not slice_entry_pattern.fullmatch(entry) for entry in entries
+ ):
+ plopm_error(
+ f"expected three i,j,k entries separated by commas, using "
+ f"positive indices, ':', or ranges, not "
+ f"{cli_error_value(f'-s {selection}')}."
+ )
+ if all(not entry for entry in entries):
+ plopm_error(
+ f"at least one slice entry must be provided with "
+ f"{cli_error_value(f'-s {selection}')}."
+ )
+ colon_entries = 0
+ for entry in entries:
+ if ":" not in entry:
+ continue
+ colon_entries += 1
+ if entry != ":":
+ start, end = (int(index) for index in entry.split(":"))
+ if start > end:
+ plopm_error(
+ f"the end of range {cli_error_value(entry)} in "
+ f"{cli_error_value(f'-s {selection}')} must not be smaller "
+ "than the start."
+ )
+ if colon_entries > 1:
+ plopm_error(
+ f"only one slice direction in "
+ f"{cli_error_value(f'-s {selection}')} can contain ':' or an index "
+ "range."
+ )
+ slice_entries.append(entries)
+
+ restart = cmdargs.restart
+ restart_pattern = re.compile(
+ rf"(?:-1|"
+ rf"{non_negative_integer}(?:,{non_negative_integer})*|"
+ rf"{non_negative_integer}:{non_negative_integer}"
+ rf"(?::{positive_integer})?)"
+ )
+ if not restart_pattern.fullmatch(restart):
+ plopm_error(
+ f"expected '-1', non-negative restart indices separated by "
+ f"commas, or 'start:end[:step]', not "
+ f"{cli_error_value(f'-r {restart}')}."
+ )
+
+ if ":" in restart:
+ restart_range = [int(value) for value in restart.split(":")]
+ if restart_range[0] > restart_range[1]:
+ plopm_error(
+ f"the end of restart range {cli_error_value(f'-r {restart}')} must "
+ "not be smaller than the start."
+ )
+
+ list_options = [
+ ("-c", "colors"),
+ ("-ls", "linestyle"),
+ ]
+ for option, name in list_options:
+ value = getattr(cmdargs, name)
+ if value and any(not entry for entry in value.split(",")):
+ plopm_error(
+ f"entries in {cli_error_value(f'{option} {value}')} cannot be empty."
+ )
+
+ line_widths = cmdargs.linewidth
+ if line_widths:
+ width_values = __parse_number_list("-lw", line_widths)
+ if any(width <= 0 for width in width_values):
+ plopm_error(
+ f"line widths must be positive, not "
+ f"{cli_error_value(f'-lw {line_widths}')}."
+ )
+
+ remove = cmdargs.hide_map_elements
+ remove_entries = remove.split(",")
+ if len(remove_entries) != 4 or any(
+ entry not in ["0", "1"] for entry in remove_entries
+ ):
+ plopm_error(
+ f"expected four values containing only 0 or 1, not "
+ f"{cli_error_value(f'-hide {remove}')}."
+ )
+
+ subfigs = cmdargs.subplot_grid
+ if subfigs:
+ subfig_entries = subfigs.split(",")
+ if len(subfig_entries) != 2 or any(
+ not re.fullmatch(positive_integer, entry) for entry in subfig_entries
+ ):
+ plopm_error(
+ f"expected two positive integers separated by a comma, such "
+ f"as {cli_error_value('-sg 2,2')}, not "
+ f"{cli_error_value(f'-sg {subfigs}')}."
+ )
+
+ colorbar_axis = cmdargs.colorbar_position
+ if colorbar_axis != "empty":
+ colorbar_axis_values = __parse_number_list(
+ "-cbp",
+ colorbar_axis,
+ 4,
+ )
+ if colorbar_axis_values[0] < 0 or colorbar_axis_values[1] < 0:
+ plopm_error(
+ f"the left and bottom positions in "
+ f"{cli_error_value(f'-cbp {colorbar_axis}')} cannot be negative."
+ )
+ if colorbar_axis_values[2] <= 0 or colorbar_axis_values[3] <= 0:
+ plopm_error(
+ f"the width and height in "
+ f"{cli_error_value(f'-cbp {colorbar_axis}')} must be positive."
+ )
+
+ grid = cmdargs.grid_edges
+ if grid:
+ grid_entries = grid.split(",")
+ if len(grid_entries) != 2 or not grid_entries[0] or not grid_entries[1]:
+ plopm_error(
+ f"expected a color and line width separated by a comma, not "
+ f"{cli_error_value(f'-ge {grid}')}."
+ )
+ if _parse_number("-ge", grid_entries[1]) < 0:
+ plopm_error(
+ f"the line width in {cli_error_value(f'-ge {grid}')} cannot be "
+ "negative."
+ )
+
+ csv_columns = cmdargs.csv_columns
+ if csv_columns:
+ csv_specifications = csv_columns.split(";")
+ for specification in csv_specifications:
+ if not specification:
+ continue
+ column_entries = specification.split(",")
+ if len(column_entries) not in [2, 3] or any(
+ not re.fullmatch(positive_integer, entry) for entry in column_entries
+ ):
+ plopm_error(
+ f"each non-empty specification in "
+ f"{cli_error_value(f'-cc {csv_columns}')} must contain two "
+ "column indices for a time series or three column indices "
+ "for a spatial map."
+ )
+ if len(set(column_entries)) != len(column_entries):
+ plopm_error(
+ f"column indices within each specification in "
+ f"{cli_error_value(f'-cc {csv_columns}')} must be different."
+ )
+
+ histogram = cmdargs.histogram
+ if histogram:
+ histogram_specifications = histogram.split()
+ for specification in histogram_specifications:
+ histogram_entries = specification.split(",")
+ if len(histogram_entries) not in [1, 2]:
+ plopm_error(
+ f"expected 'bins', 'bins,norm', or 'bins,lognorm', not "
+ f"{cli_error_value(f'-hist {specification}')}."
+ )
+ if not re.fullmatch(
+ positive_integer,
+ histogram_entries[0],
+ ):
+ plopm_error(
+ f"the number of bins in "
+ f"{cli_error_value(f'-hist {specification}')} must be a positive "
+ "integer."
+ )
+ if len(histogram_entries) == 2 and histogram_entries[1] not in [
+ "norm",
+ "lognorm",
+ ]:
+ plopm_error(
+ f"the distribution in "
+ f"{cli_error_value(f'-hist {specification}')} must be 'norm' or "
+ "'lognorm'."
+ )
+
+ band_properties = cmdargs.fill_between_style
+ if band_properties:
+ band_entries = band_properties.split(",")
+ if len(band_entries) % 2 != 0 or any(not color for color in band_entries[::2]):
+ plopm_error(
+ f"expected color and alpha pairs, not "
+ f"{cli_error_value(f'-fb {band_properties}')}."
+ )
+ try:
+ alpha_values = [float(alpha) for alpha in band_entries[1::2]]
+ except ValueError:
+ alpha_values = []
+ if not alpha_values or any(alpha < 0 or alpha > 1 for alpha in alpha_values):
+ plopm_error(
+ f"alpha values in {cli_error_value(f'-fb {band_properties}')} must "
+ "be between 0 and 1."
+ )
+ if cmdargs.ensemble not in ["1", "3"]:
+ plopm_error(
+ f"{cli_error_value('-fb')} can only be used with "
+ f"{cli_error_value('-ens 1')} or {cli_error_value('-ens 3')}."
+ )
+
+ log_values = cmdargs.color_log.split(",")
+
+ if cmdargs.color_log_ticks and "1" not in log_values:
+ plopm_error(
+ f"{cli_error_value('-clt')} requires at least one logarithmic color "
+ f"scale enabled with {cli_error_value('-clog')}."
+ )
+
+ if cmdargs.mask_threshold != "1e-3" and not cmdargs.mask_variable:
+ plopm_error(
+ f"{cli_error_value('-mt')} can only be changed when "
+ f"{cli_error_value('-mv')} is used."
+ )
+
+ if (
+ cmdargs.distance
+ and "sensor" in cmdargs.distance
+ and any(
+ any(not re.fullmatch(positive_integer, entry) for entry in entries)
+ for entries in slice_entries
+ )
+ ):
+ plopm_error(
+ f"a sensor distance requires each location provided with "
+ f"{cli_error_value('-s')} to contain three positive indices."
+ )
+
+ vtk_names = cmdargs.vtk_names
+ if vtk_names:
+ vtk_name_entries = vtk_names.split(",")
+ if any(not name for name in vtk_name_entries):
+ plopm_error(
+ f"VTK variable names in {cli_error_value(f'-vn {vtk_names}')} "
+ "cannot be empty."
+ )
+
+ valid_vtk_formats = [
+ "Float64",
+ "Float32",
+ "Float16",
+ "Int64",
+ "UInt64",
+ "Int32",
+ "UInt32",
+ "Int16",
+ "UInt16",
+ "Int8",
+ "UInt8",
+ ]
+ vtk_formats = cmdargs.vtk_format.split(",")
+ if any(vtk_format not in valid_vtk_formats for vtk_format in vtk_formats):
+ plopm_error(
+ f"expected VTK formats from {', '.join(valid_vtk_formats)}, not "
+ f"{cli_error_value(f'-vf {cmdargs.vtk_format}')}."
+ )
+
+ vtk_options = {
+ "-fp": ("flow_path", "flow"),
+ "-vf": ("vtk_format", "Float64"),
+ "-vn": ("vtk_names", ""),
+ }
+ if not vtk_mode:
+ invalid_options = [
+ option
+ for option, (name, default) in vtk_options.items()
+ if getattr(cmdargs, name) != default
+ ]
+ if invalid_options:
+ formatted_options = ", ".join(
+ cli_error_value(option) for option in invalid_options
+ )
+ plopm_error(
+ f"{formatted_options} can only be used with "
+ f"{cli_error_value('-m vtk')}, not {cli_error_value(f'-m {mode}')}."
+ )
+ else:
+ try:
+ flow_arguments = shlex.split(cmdargs.flow_path)
+ except ValueError:
+ flow_arguments = []
+
+ if not flow_arguments:
+ plopm_error(
+ f"the OPM Flow command "
+ f"{cli_error_value(f'-fp {cmdargs.flow_path}')} cannot be empty."
+ )
+
+ try:
+ flow_result = subprocess.run(
+ [*flow_arguments, "-h"],
+ stdout=subprocess.DEVNULL,
+ stderr=subprocess.STDOUT,
+ check=False,
+ )
+ except OSError:
+ flow_result = None
+
+ if flow_result is None or flow_result.returncode != 0:
+ plopm_error(
+ f"the OPM Flow executable "
+ f"{cli_error_value(f'-fp {cmdargs.flow_path}')} is not available or "
+ "not working."
+ )
+
+ if not gif_mode:
+ gif_options = {
+ "-gi": ("gif_interval", "1000"),
+ "-gl": ("gif_loop", "0"),
+ }
+ invalid_options = [
+ option
+ for option, (name, default) in gif_options.items()
+ if getattr(cmdargs, name) != default
+ ]
+ if invalid_options:
+ formatted_options = ", ".join(
+ cli_error_value(option) for option in invalid_options
+ )
+ plopm_error(
+ f"{formatted_options} can only be used with "
+ f"{cli_error_value('-m gif')}, not {cli_error_value(f'-m {mode}')}."
+ )
+
+
+def _parse_number(option: str, value: str) -> float:
+ """Parse one numeric command-line value.
+
+ Parameters
+ ----------
+ option : str
+ Option name used in an error message.
+ value : str
+ Value to convert.
+
+ Returns
+ -------
+ float
+ Parsed numeric value.
+
+ Raises
+ ------
+ SystemExit
+ If the value is not numeric.
+
+ """
+ try:
+ number = float(value)
+ except ValueError:
+ plopm_error(f"expected a number, not {cli_error_value(f'{option} {value}')}.")
+ return number
+
+
+def __parse_number_list(
+ option: str,
+ value: str,
+ expected_length: int | None = None,
+) -> list:
+ """Parse comma-separated numeric values.
+
+ Parameters
+ ----------
+ option : str
+ Option name used in an error message.
+ value : str
+ Comma-separated values to convert.
+ expected_length : int, optional
+ Required number of values.
+
+ Returns
+ -------
+ list[float]
+ Parsed numeric values.
+
+ Raises
+ ------
+ SystemExit
+ If a value is not numeric or the length is invalid.
+
+ """
+ entries = value.split(",")
+ if expected_length is not None and len(entries) != expected_length:
+ plopm_error(
+ f"expected {expected_length} numbers separated by commas, "
+ f"not {cli_error_value(f'{option} {value}')}."
+ )
+ try:
+ numbers = [float(entry) for entry in entries]
+ except ValueError:
+ plopm_error(
+ "expected numbers separated by commas, "
+ f"not {cli_error_value(f'{option} {value}')}."
+ )
+ return numbers
+
+
+
+
+
+# SPDX-FileCopyrightText: 2024-2026 NORCE Research AS
+# SPDX-License-Identifier: GPL-3.0
+# pylint: disable=W0123,R0915,R0912,R1702,R0914,R0916
+
+"""Build and normalize configuration for plopm workflows.
+
+The module converts parsed CLI arguments into :class:`PlopmConfig`, discovers
+simulation cases, normalizes per-variable plotting settings, selects summary or
+spatial processing, and defines unit conversions used by the readers.
+"""
+
+import argparse
+import copy
+import os
+import shutil
+import sys
+from typing import cast
+
+import matplotlib
+import matplotlib.pyplot as plt
+from opm.io.ecl import EclFile as OpmFile
+from opm.io.ecl import ESmry as OpmSummary
+
+from plopm.config.config import PlopmConfig
+from plopm.utils.terminal import (
+ cli_current_value,
+ cli_error_value,
+ cli_info_value,
+ plopm_error,
+ plopm_info,
+)
+
+
+
+[docs]
+def build_config(cmdargs: argparse.Namespace) -> PlopmConfig:
+ """Build a plopm configuration from parsed CLI arguments.
+
+ The function expands case and difference-input paths, parses list-like
+ options, normalizes slice and restart selections, and initializes plotting
+ defaults shared by summary, map, and VTK workflows.
+
+ Parameters
+ ----------
+ cmdargs : argparse.Namespace
+ Command-line arguments returned by the plopm parser.
+
+ Returns
+ -------
+ PlopmConfig
+ Parsed and partially normalized runtime configuration.
+
+ """
+
+ cfg = PlopmConfig()
+ cfg.output_dir = os.path.abspath(cmdargs.output_dir)
+ names = cmdargs.input.split(" ")
+ names = [var.split(" ") for var in names]
+ cfg.case_labels = names
+
+ for name in ["gif", "csv", "png", "vtk"]:
+ setattr(cfg, name, cmdargs.format == name)
+
+ cfg.difference_input = cmdargs.difference_input
+ cfg.ensemble = int(cmdargs.ensemble)
+
+ if cfg.difference_input:
+ if cfg.difference_input[-1] in [".", "/"]:
+ cfg.difference_input = _find_first_case(cfg.difference_input, ".EGRID")
+ if names[0][0][-1] in [".", "/"]:
+ names[0][0] = _find_first_case(names[0][0], ".EGRID")
+ elif names[0][0][-1] in [".", "/"]:
+ folders = names[0]
+ names = []
+ for index, folder in enumerate(folders):
+ if cfg.ensemble > 0 or index == 0:
+ names.append([])
+ if cfg.vtk:
+ names[-1] = _find_all_cases(folder, ".DATA")
+ else:
+ names[-1] = _find_all_cases(folder, ".SMSPEC")
+
+ cfg.cases = names
+ cfg.case = names[0][0]
+ cfg.variables = cmdargs.variable.lower().split(",")
+ _join_block_vars(cfg)
+ cfg.stress_coefficient = float(cmdargs.stress_coefficient)
+
+ for cfg_name, cmdarg_name in [
+ ("vtk_names", "vtk_names"),
+ ("filename", "filename"),
+ ]:
+ setattr(cfg, cfg_name, getattr(cmdargs, cmdarg_name).split(" "))
+
+ cfg.mass_vars = ["gasm", "dism", "liqm", "vapm", "co2m", "h2om"]
+ cfg.mass_fracs = ["xco2l", "xh2ov", "xco2v", "xh2ol"]
+ cfg.caprock_vars = ["limipres", "overpres", "objepres"]
+ for cfg_name, cmdarg_name in [
+ ("filters", "filters"),
+ ("restart", "restart"),
+ ("scale_factor", "scale_factor"),
+ ("vtk_format", "vtk_format"),
+ ]:
+ setattr(cfg, cfg_name, getattr(cmdargs, cmdarg_name).split(","))
+
+ if cfg.restart[0] == "-1":
+ cfg.restart = [-1]
+ elif ":" in cfg.restart[0]:
+ cfg.rst_range = True
+ vals = cfg.restart[0].split(":")
+ if len(vals) == 3:
+ cfg.restart = list(
+ range(
+ int(vals[0]),
+ int(vals[1]) + 1,
+ int(vals[2]),
+ )
+ )
+ else:
+ cfg.restart = list(
+ range(
+ int(vals[0]),
+ int(vals[1]) + 1,
+ )
+ )
+ if cfg.filename[0]:
+ width = len(str(cfg.restart[-1]))
+ cfg.filename = [
+ cfg.filename[0] + f"{restart_value}".zfill(width)
+ for restart_value in cfg.restart
+ ]
+ else:
+ if "," in cmdargs.restart and (cfg.png or cfg.csv):
+ cfg.rst_range = True
+ width = len(str(cfg.restart[-1]))
+ cfg.filename = [
+ cfg.filename[0] + f"{restart_value}".zfill(width)
+ for restart_value in cfg.restart
+ ]
+ cfg.restart = [int(restart_value) for restart_value in cfg.restart]
+ for name in ["vtk_format", "scale_factor", "vtk_names"]:
+ if len(getattr(cfg, name)) < len(cfg.variables):
+ setattr(
+ cfg,
+ name,
+ [getattr(cfg, name)[0]] * len(cfg.variables),
+ )
+ if not os.path.exists(cfg.output_dir):
+ os.makedirs(cfg.output_dir, exist_ok=True)
+ if cfg.vtk:
+ return cfg
+
+ cfg.csv_columns = cmdargs.csv_columns.split(";")
+ cfg.csv_columns = [
+ [int(val) if val else "" for val in var.split(",")] for var in cfg.csv_columns
+ ]
+
+ allcsvs = True
+ for val in cfg.csv_columns:
+ if not val[0]:
+ allcsvs = False
+ elif len(val) == 2:
+ cfg.csv_column_summary = True
+
+ if allcsvs:
+ cfg.variables = ["csv"]
+
+ max_count = max(len(cfg.cases[0]), len(cfg.variables))
+ if len(cfg.csv_columns) == 1 and not cfg.csv_columns[0][0]:
+ cfg.csv_columns = [cfg.csv_columns[0]] * (max_count + 1)
+
+ for cfg_name, cmdarg_name in [
+ ("mask_variable", "mask_variable"),
+ ("linewidth", "linewidth"),
+ ("linestyle", "linestyle"),
+ ("inactive_color", "inactive_color"),
+ ]:
+ setattr(cfg, cfg_name, getattr(cmdargs, cmdarg_name).lower())
+
+ for cfg_name, cmdarg_name in [
+ ("fontsize", "fontsize"),
+ ("mask_threshold", "mask_threshold"),
+ ("gif_interval", "gif_interval"),
+ ]:
+ setattr(cfg, cfg_name, float(getattr(cmdargs, cmdarg_name)))
+
+ for cfg_name, cmdarg_name in [
+ ("colorbar_ticks", "colorbar_ticks"),
+ ("title", "title"),
+ ]:
+ setattr(cfg, cfg_name, getattr(cmdargs, cmdarg_name).split(" "))
+
+ for cfg_name, cmdarg_name in [
+ ("clim", "clim"),
+ ("translation", "translation"),
+ ("histogram", "histogram"),
+ ]:
+ setattr(cfg, cfg_name, getattr(cmdargs, cmdarg_name).split(" "))
+
+ for cfg_name, cmdarg_name in [
+ ("suptitle", "suptitle"),
+ ("fill_between_style", "fill_between_style"),
+ ("colorbar_label", "colorbar_label"),
+ ]:
+ setattr(cfg, cfg_name, getattr(cmdargs, cmdarg_name))
+
+ cfg.clim = [var.split(",") for var in cfg.clim]
+ cfg.translation = [var.split(",") for var in cfg.translation]
+ cfg.colors_raw = cmdargs.colors
+ cfg.colorbar_format = cmdargs.colorbar_format
+ cfg.fc = cmdargs.facecolor
+ cfg.legend_labels = cmdargs.legend_labels.split(" ")
+ cfg.legend_labels = [var.split(" ") for var in cfg.legend_labels]
+ cfg.hide_map_elements = [int(val) for val in cmdargs.hide_map_elements.split(",")]
+ cfg.global_range = int(cmdargs.global_range) == 1
+
+ for cfg_name, cmdarg_name in [
+ ("equal_aspect", "equal_aspect"),
+ ("remove_duplicate_labels", "remove_duplicate_labels"),
+ ("gif_loop", "gif_loop"),
+ ("list_variables", "list_variables"),
+ ("step_plot", "step_plot"),
+ ]:
+ setattr(cfg, cfg_name, int(getattr(cmdargs, cmdarg_name)) == 1)
+
+ for cfg_name, cmdarg_name in [
+ ("figsize", "figsize"),
+ ("distance", "distance"),
+ ("aggregation", "aggregation"),
+ ("rotation", "rotation"),
+ ("color_log", "color_log"),
+ ("legend_location", "legend_location"),
+ ("axis_grid", "axis_grid"),
+ ]:
+ setattr(cfg, cfg_name, getattr(cmdargs, cmdarg_name).split(","))
+
+ for cfg_name, cmdarg_name in [
+ ("dpi", "dpi"),
+ ("time_units", "time_units"),
+ ("colorbar_tick_count", "colorbar_tick_count"),
+ ("grid_edges", "grid_edges"),
+ ]:
+ setattr(cfg, cfg_name, getattr(cmdargs, cmdarg_name).split(","))
+
+ for cfg_name, cmdarg_name in [
+ ("dual_grid", "dual_grid"),
+ ("subplot_grid", "subplot_grid"),
+ ("min_threshold", "min_threshold"),
+ ("max_threshold", "max_threshold"),
+ ]:
+ setattr(cfg, cfg_name, getattr(cmdargs, cmdarg_name).split(","))
+
+ for axis_name in ["x", "y"]:
+ setattr(
+ cfg,
+ f"{axis_name}units",
+ getattr(cmdargs, f"{axis_name}units"),
+ )
+ setattr(
+ cfg,
+ f"{axis_name}label",
+ getattr(cmdargs, f"{axis_name}label").split(" "),
+ )
+ setattr(
+ cfg,
+ f"{axis_name}format",
+ getattr(cmdargs, f"{axis_name}format").split(","),
+ )
+ setattr(
+ cfg,
+ f"{axis_name}tick_count",
+ getattr(cmdargs, f"{axis_name}tick_count").split(","),
+ )
+ setattr(
+ cfg,
+ f"{axis_name}log",
+ getattr(cmdargs, f"{axis_name}log").split(","),
+ )
+ setattr(
+ cfg,
+ f"{axis_name}lim",
+ getattr(cmdargs, f"{axis_name}lim").split(" "),
+ )
+ setattr(
+ cfg,
+ f"{axis_name}lim",
+ [var.split(",") for var in getattr(cfg, f"{axis_name}lim")],
+ )
+
+ if cmdargs.color_log_ticks:
+ cfg.color_log_ticks = [
+ float(val) for val in cmdargs.color_log_ticks[1:-1].split(",")
+ ]
+
+ if cfg.colorbar_ticks[0]:
+ for index, values in enumerate(cfg.colorbar_ticks):
+ cfg.colorbar_ticks[index] = [val.strip() for val in values[1:-1].split(",")]
+ if cmdargs.colorbar_position != "empty":
+ cfg.colorbar_position = cast(
+ tuple[float, float, float, float],
+ tuple(map(float, cmdargs.colorbar_position.split(","))),
+ )
+
+ cfg.slice = cmdargs.slice.split(" ")
+ cfg.slice = [
+ [val if val else [-2, -2] for val in var.split(",")] for var in cfg.slice
+ ]
+ if [-2, -2] in cfg.slice[0]:
+ for slice_index, var in enumerate(cfg.slice):
+ for value_index, val in enumerate(var):
+ if val[0] != -2:
+ if val == ":":
+ pass
+ elif ":" in val:
+ vals = val.split(":")
+ cfg.slice[slice_index][value_index] = [
+ int(vals[0]) - 1,
+ int(vals[1]),
+ ]
+ else:
+ int_value = int(val)
+ cfg.slice[slice_index][value_index] = [int_value - 1, int_value]
+ elif ":" in cfg.slice[0]:
+ cfg.layer = True
+ for slice_index, var in enumerate(cfg.slice):
+ for value_index, val in enumerate(var):
+ if val != ":":
+ cfg.slice[slice_index][value_index] = int(val) - 1
+ else:
+ cfg.slice[slice_index][value_index] = -1
+ else:
+ cfg.sensor = True
+ for slice_index, var in enumerate(cfg.slice):
+ for value_index, val in enumerate(var):
+ cfg.slice[slice_index][value_index] = int(val) - 1
+
+ cfg.summary_mass = ["fwcdm", "fgipm"]
+
+ cfg.colors_default = [
+ "k",
+ "b",
+ "#ff7f0e",
+ "#2ca02c",
+ "#d62728",
+ "#9467bd",
+ "#8c564b",
+ "#e377c2",
+ "#7f7f7f",
+ "#bcbd22",
+ "#17becf",
+ "#1f77b4",
+ "r",
+ ]
+
+ cfg.linestyle_default = [
+ "-",
+ "--",
+ (0, (1, 1)),
+ "-.",
+ (0, (1, 10)),
+ (0, (1, 1)),
+ (5, (10, 3)),
+ (0, (5, 10)),
+ (0, (5, 5)),
+ (0, (5, 1)),
+ (0, (3, 10, 1, 10)),
+ (0, (3, 5, 1, 5)),
+ (0, (3, 1, 1, 1)),
+ (0, (3, 5, 1, 5, 1, 5)),
+ (0, (3, 10, 1, 10, 1, 10)),
+ (0, (3, 1, 1, 1, 1, 1)),
+ (0, ()),
+ ]
+ for val in cfg.variables:
+ for oper in ["=", "<", ">"]:
+ if oper in val:
+ cfg.discrete = False
+
+ cfg.linewidth_values = ["1"] * len(cfg.cases[0])
+
+ font = {"family": "normal", "weight": "normal", "size": cfg.fontsize}
+ matplotlib.rc("font", **font)
+ plt.rcParams.update(
+ {
+ "text.usetex": shutil.which("latex") is not None,
+ "font.family": "monospace",
+ "legend.columnspacing": 0.9,
+ "legend.handlelength": 3.5,
+ "legend.fontsize": cfg.fontsize,
+ "lines.linewidth": 4,
+ "axes.titlesize": cfg.fontsize,
+ "axes.grid": False,
+ "figure.figsize": (float(cfg.figsize[0]), float(cfg.figsize[1])),
+ }
+ )
+
+ if len(cfg.filename) < len(cfg.variables):
+ cfg.filename = [cfg.filename[0]] * len(cfg.variables)
+
+ if len(cfg.clim) < len(cfg.variables):
+ cfg.clim = [cfg.clim[0]] * len(cfg.variables)
+
+ if cfg.difference_input and len(cfg.rotation) < 2:
+ cfg.rotation = [cfg.rotation[0]] * 2
+ elif len(cfg.rotation) < len(cfg.cases[0]):
+ cfg.rotation = [cfg.rotation[0]] * len(cfg.cases[0])
+
+ if cfg.difference_input and len(cfg.translation) < 2:
+ cfg.translation = [cfg.translation[0]] * 2
+
+ if len(cfg.translation) < len(cfg.cases[0]):
+ cfg.translation = [cfg.translation[0]] * len(cfg.cases[0])
+
+ if cfg.difference_input and len(cfg.slice) < 2:
+ cfg.slice = [cfg.slice[0]] * 2
+
+ for val in [
+ "aggregation",
+ "filters",
+ "colorbar_ticks",
+ "csv_columns",
+ "dual_grid",
+ "slice",
+ "title",
+ ]:
+ if len(getattr(cfg, val)) < max_count:
+ if val == "slice":
+ current = getattr(cfg, val)
+ setattr(
+ cfg,
+ val,
+ [copy.deepcopy(current[0]) for _ in range(max_count)],
+ )
+ else:
+ setattr(cfg, val, [getattr(cfg, val)[0]] * max_count)
+ elif len(cfg.restart) > 1 and cfg.subplot_grid[0]:
+ if (
+ len(getattr(cfg, val)) >= max(max_count, len(cfg.restart))
+ and val == "title"
+ ):
+ continue
+ if val == "slice":
+ if cfg.gif and len(cfg.slice) >= len(cfg.cases[0]):
+ continue
+ current = getattr(cfg, val)
+ setattr(
+ cfg,
+ val,
+ [copy.deepcopy(current[0]) for _ in range(len(cfg.restart))],
+ )
+ else:
+ setattr(cfg, val, [getattr(cfg, val)[0]] * len(cfg.restart))
+
+ if len(cfg.restart) > 1 and cfg.subplot_grid[0]:
+ cfg.filename = [cmdargs.filename]
+ if cfg.difference_input:
+ cfg.aggregation = [cfg.aggregation[0]] * 2
+ cfg.filters = [cfg.filters[0]] * 2
+
+ for val in [
+ "xformat",
+ "yformat",
+ "xlog",
+ "ylog",
+ "xlabel",
+ "ylabel",
+ "legend_labels",
+ "time_units",
+ "legend_location",
+ "dpi",
+ "ytick_count",
+ "xtick_count",
+ "filename",
+ "axis_grid",
+ "colorbar_tick_count",
+ "color_log",
+ "min_threshold",
+ "max_threshold",
+ ]:
+ if len(getattr(cfg, val)) < len(cfg.variables):
+ setattr(cfg, val, [getattr(cfg, val)[0]] * len(cfg.variables))
+
+ return cfg
+
+
+
+def _find_all_cases(folder: str, suffix: str) -> list:
+ """Find all simulation cases below a folder.
+
+ Parameters
+ ----------
+ folder : str
+ Folder to search recursively.
+ suffix : str
+ File suffix identifying a simulation case.
+
+ Returns
+ -------
+ list[str]
+ Sorted case paths without the identifying suffix.
+
+ """
+ folder_path = folder
+ if folder_path[0] != ".":
+ folder_path = "./" + folder_path
+ cases_found = []
+ for root, _, files in os.walk(folder_path):
+ for file in files:
+ if file.endswith(suffix):
+ cases_found.append(os.path.join(root, file)[2 : -len(suffix)])
+ return sorted(cases_found)
+
+
+def _find_first_case(folder: str, suffix: str) -> str:
+ """Find the first simulation case below a folder.
+
+ Parameters
+ ----------
+ folder : str
+ Folder to search recursively.
+ suffix : str
+ File suffix identifying a simulation case.
+
+ Returns
+ -------
+ str
+ First case path without the suffix, or the input folder when no
+ matching file is found.
+
+ """
+ folder_path = folder
+ if folder_path[0] != ".":
+ folder_path = "./" + folder_path
+ for root, _, files in os.walk(folder_path):
+ for file in files:
+ if file.endswith(suffix):
+ return os.path.join(root, file)[2 : -len(suffix)]
+ return folder
+
+
+def _join_block_vars(cfg: PlopmConfig) -> None:
+ """Rejoin comma-separated indices in block variables.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Configuration whose variable expressions are updated in place.
+
+ """
+ vrs_in = cfg.variables
+ count = len(vrs_in)
+ variables = []
+ index = 0
+ while index < count:
+ if index < count - 2 and ":" in vrs_in[index] and vrs_in[index + 1].isnumeric():
+ variables.append(
+ vrs_in[index] + "," + vrs_in[index + 1] + "," + vrs_in[index + 2]
+ )
+ index += 3
+ else:
+ variables.append(vrs_in[index])
+ index += 1
+ cfg.variables = variables
+
+
+
+[docs]
+def init_maps(cfg: PlopmConfig) -> None:
+ """Normalize settings used by spatial maps.
+
+ The function selects default units, colorbar formats, and colormaps; expands
+ per-variable limits and formats; and initializes spatial coordinate scales.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Configuration updated in place for map generation.
+
+ """
+ cfg.units = [" [-]", " [mD]", " [mD]", r" [m$^3$]", " [-]", " [-]"]
+ cfg.cb_formats = [".1f", ".0f", ".0f", ".2e", ".0f", ".0f"]
+ cfg.colormaps = ["jet", "turbo", "turbo", "terrain", "tab20b", "tab20b"]
+ cmdisc = [
+ "Pastel1",
+ "Pastel2",
+ "Paired",
+ "Accent",
+ "Dark2",
+ "Set1",
+ "Set2",
+ "Set3",
+ "tab10",
+ "tab20",
+ "tab20b",
+ "tab20c",
+ "cet_glasbey_bw",
+ "cet_glasbey",
+ "cet_glasbey_cool",
+ "cet_glasbey_warm",
+ "cet_glasbey_dark",
+ "cet_glasbey_light",
+ "cet_glasbey_category10",
+ "cet_glasbey_hv",
+ ]
+ cfg.disc_colormaps = [cmap + "_r" for cmap in cmdisc] + cmdisc
+ if cfg.colors_raw:
+ cfg.colormaps = cfg.colors_raw.split(",")
+ elif cfg.difference_input:
+ cfg.colormaps = ["RdYlGn"]
+ elif cfg.mask_variable:
+ cfg.colormaps = ["RdGy_r"]
+ variables = cfg.variables
+ if variables:
+ first_variable = variables[0]
+ if first_variable in ["wells", "faults"]:
+ if cfg.aggregation[0]:
+ if cfg.aggregation[0] not in ["min", "max"]:
+ plopm_error(
+ f"Unsuported value {cli_error_value(f'-agg {cfg.aggregation[0]}')} for "
+ f"{cli_info_value(f'-v {first_variable}')}. Supported values are "
+ f"{cli_current_value('-agg min')} and {cli_current_value('-agg max')}."
+ )
+ cfg.slice_mode = cfg.aggregation[0]
+ else:
+ cfg.slice_mode = "min"
+ if not cfg.colors_raw:
+ cfg.units = [" [-]"]
+ cfg.colormaps = ["nipy_spectral"]
+ cfg.cb_formats = [".0f"]
+ if (
+ "num" in first_variable
+ and not cfg.mask_variable
+ and not cfg.difference_input
+ and not cfg.colors_raw
+ ):
+ cfg.colormaps = ["tab20"]
+ cfg.units = [" [-]"]
+ cfg.cb_formats = [".0f"]
+ if "index" in first_variable:
+ cfg.units = [" [-]"]
+ cfg.cb_formats = [".0f"]
+ if cfg.colorbar_format:
+ cfg.cb_formats = cfg.colorbar_format.split(",")
+ elif len(variables) == 1 and "num" in variables[0]:
+ cfg.cb_formats = [".0f"]
+ elif cfg.difference_input:
+ cfg.cb_formats = [".1e"]
+ nvars = len(variables)
+ if len(cfg.colormaps) < nvars or (
+ nvars == len(cfg.cases[0]) and len(cfg.cases[0]) > 1 and not cfg.colors_raw
+ ):
+ cfg.colormaps = [cfg.colormaps[0]] * nvars
+ if len(cfg.xlim) < nvars:
+ cfg.xlim = [cfg.xlim[0]] * nvars
+ if len(cfg.ylim) < nvars:
+ cfg.ylim = [cfg.ylim[0]] * nvars
+ if len(cfg.cb_formats) < nvars:
+ cfg.cb_formats = [cfg.cb_formats[0]] * nvars
+ cfg.xscale, cfg.xunit = spatial_unit(cfg.xunits)
+ cfg.yscale, cfg.yunit = spatial_unit(cfg.yunits)
+
+
+
+
+[docs]
+def spatial_unit(unit: str) -> tuple[float, str]:
+ """Get the conversion and label for a spatial unit.
+
+ Parameters
+ ----------
+ unit : str
+ Spatial-unit code.
+
+ Returns
+ -------
+ tuple[float, str]
+ Factor converting metres and the formatted unit label.
+
+ """
+ return {
+ "m": (1.0, " [m]"),
+ "km": (1e-3, " [km]"),
+ "cm": (1e2, " [cm]"),
+ "mm": (1e3, " [mm]"),
+ }.get(unit, (1.0, ""))
+
+
+
+
+[docs]
+def mass_unit(mskl: float) -> str:
+ """Get the display unit for a mass scale factor.
+
+ Parameters
+ ----------
+ mskl : float
+ Factor applied to quantities stored in kilograms.
+
+ Returns
+ -------
+ str
+ Matplotlib-formatted mass unit, or an empty string when unknown.
+
+ """
+ return {
+ 1e-3: " [t]",
+ 1e-6: " [Kt]",
+ 1e-9: " [Mt]",
+ 1e3: " [g]",
+ 1e6: " [mg]",
+ 1: " [kg]",
+ }.get(mskl, "")
+
+
+
+
+[docs]
+def is_summary(cfg: PlopmConfig) -> bool:
+ """Determine whether the request uses one-dimensional output.
+
+ The decision considers explicit series options, special tabulated
+ properties, summary-vector availability, and requests to list variables.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized configuration and primary case path.
+
+ Returns
+ -------
+ bool
+ ``True`` when the request should use the summary plotting workflow.
+
+ """
+ name = cfg.case
+ variables = cfg.variables
+ first_variable = variables[0] if variables else ""
+ ntot = 0
+ if cfg.list_variables:
+ for ext in ["INIT", "UNRST"]:
+ file = f"{name}.{ext}"
+ if os.path.isfile(file):
+ reader = OpmFile(file)
+ keys = [
+ var[0]
+ for var in reader.arrays
+ if var[0]
+ not in ["INTEHEAD", "LOGIHEAD", "DOUBHEAD", "TABDIMS", "TAB"]
+ ]
+ if ext == "UNRST":
+ ntot = reader.count("PRESSURE")
+ plopm_info(
+ f"the available {cli_info_value('-v')} variables for "
+ f"{cli_info_value(file)} are:"
+ )
+ print(keys)
+ if ext == "UNRST":
+ plopm_info(
+ f"the available {cli_info_value('-r')} restarts for "
+ f"{cli_info_value(file)} are:"
+ )
+ print(list(range(ntot)))
+ if (
+ cfg.sensor
+ or cfg.layer
+ or cfg.distance[0]
+ or cfg.histogram[0]
+ or cfg.csv_column_summary
+ ):
+ return True
+ if (
+ first_variable[:3] in ["krw", "krg"]
+ or first_variable[:4] in ["krow", "krog", "pcow", "pcog", "pcwg"]
+ or first_variable[:6] == "pcfact"
+ or (
+ first_variable[:8] == "permfact"
+ and cfg.slice == [[[-2, -2], [0, 1], [-2, -2]]]
+ )
+ ):
+ return True
+ smspec_file = f"{name}.SMSPEC"
+ if os.path.isfile(smspec_file):
+ summary = OpmSummary(smspec_file).keys()
+ if cfg.list_variables:
+ plopm_info(
+ f"the available {cli_info_value('-v')} variables for "
+ f"{cli_info_value(smspec_file)} are:"
+ )
+ print(summary)
+ sys.exit(0)
+ smass = cfg.summary_mass
+ for name_v in variables:
+ base = name_v.split(" ")[0].upper()
+ if base in summary or base.lower() in smass:
+ return True
+ if cfg.list_variables:
+ sys.exit(0)
+ return False
+
+
+
+
+[docs]
+def init_summary(cfg: PlopmConfig) -> None:
+ """Normalize settings used by one-dimensional plots.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Configuration updated in place with per-variable styles and labels.
+
+ """
+ variables = cfg.variables
+ nvars = len(variables)
+ cfg.ncolors = 1 if len(cfg.cases) < nvars else nvars
+ for val in ["colors_raw", "linestyle", "linewidth"]:
+ if getattr(cfg, val):
+ tmp = [var.split(",") for var in getattr(cfg, val).split(":")]
+ if len(tmp) < nvars:
+ tmp = [tmp[0]] * nvars
+ setattr(cfg, "colors" if val == "colors_raw" else val, tmp)
+ elif val == "colors_raw":
+ cfg.colors = [cfg.colors_default] * nvars
+ elif val == "linestyle":
+ cfg.linestyle = [cfg.linestyle_default] * nvars
+ else:
+ cfg.linewidth = [cfg.linewidth_values] * nvars
+ for axis_name in ["x", "y"]:
+ key = f"{axis_name}lim"
+ if len(getattr(cfg, key)) < nvars and getattr(cfg, key)[0]:
+ setattr(cfg, key, [getattr(cfg, key)[0]] * nvars)
+ if nvars == 1 and len(cfg.linewidth[0]) < len(cfg.cases[0]):
+ cfg.linewidth[0] = [cfg.linewidth[0][0]] * len(cfg.cases[0])
+ for val in [
+ "cases",
+ "title",
+ "xformat",
+ "yformat",
+ "xlog",
+ "ylog",
+ "xlabel",
+ "ylabel",
+ "legend_labels",
+ "time_units",
+ "legend_location",
+ "dpi",
+ "ytick_count",
+ "xtick_count",
+ "scale_factor",
+ "filename",
+ "axis_grid",
+ ]:
+ if len(getattr(cfg, val)) < nvars:
+ setattr(cfg, val, [getattr(cfg, val)[0]] * nvars)
+
+
+
+
+
+
+# SPDX-FileCopyrightText: 2024-2026 NORCE Research AS
+# SPDX-License-Identifier: GPL-3.0
+# pylint: disable=R1702,R0912,C0325,R0913,R0914,R0915,R0917
+
+"""Prepare slice geometry and map three-dimensional values to two dimensions.
+
+The module builds labels and coordinate meshes for xy, xz, and yz slices,
+applies optional rotation and translation, and aggregates active-cell values
+through the selected grid interval.
+"""
+
+import numpy as np
+from numpy.typing import NDArray
+
+from plopm.config.config import PlopmConfig, SimData
+from plopm.utils.readers import get_xy_coords, get_xz_coords, get_yz_coords
+
+
+
+[docs]
+def get_yz_slice(
+ cfg: PlopmConfig, data: SimData, n: int
+) -> tuple[NDArray, NDArray, str, str, int, int, str, str]:
+ """Prepare geometry and labels for a yz slice.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map and slice configuration.
+ data : SimData
+ Loaded grid data.
+ n : int
+ Slice index.
+
+ Returns
+ -------
+ tuple
+ Coordinate meshes, display and filename slice labels, mapped grid
+ dimensions, and coordinate-axis names.
+
+ """
+ slice_range = cfg.slice[n][0]
+ nx = data.nx
+ if slice_range[0] == ":":
+ cfg.slice[n][0] = [0, nx]
+ slice_title = f", slice i=0:{nx}"
+ slice_name = f"0:{nx},j,k"
+ elif slice_range[0] == slice_range[1] - 1:
+ start_index = slice_range[0] + 1
+ slice_title = f", slice i={start_index}"
+ slice_name = f"{start_index},j,k"
+ else:
+ start_index = slice_range[0] + 1
+ end_index = slice_range[1]
+ slice_title = f", slice i={start_index}:{end_index}"
+ slice_name = f"{start_index}:{end_index},j,k"
+ xc, yc = get_yz_coords(cfg, data, n)
+ mx = 2 * data.ny - 1
+ my = 2 * data.nz - 1
+ xname = "y"
+ yname = "z"
+ return xc, yc, slice_title, slice_name, mx, my, xname, yname
+
+
+
+
+[docs]
+def get_xz_slice(
+ cfg: PlopmConfig, data: SimData, n: int
+) -> tuple[NDArray, NDArray, str, str, int, int, str, str]:
+ """Prepare geometry and labels for an xz slice.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map and slice configuration.
+ data : SimData
+ Loaded grid data.
+ n : int
+ Slice index.
+
+ Returns
+ -------
+ tuple
+ Coordinate meshes, display and filename slice labels, mapped grid
+ dimensions, and coordinate-axis names.
+
+ """
+ slice_range = cfg.slice[n][1]
+ ny = data.ny
+ if slice_range[0] == ":":
+ cfg.slice[n][1] = [0, ny]
+ slice_title = f", slice j=0:{ny}"
+ slice_name = f"i,0:{ny},k"
+ elif slice_range[0] == slice_range[1] - 1:
+ start_index = slice_range[0] + 1
+ slice_title = f", slice j={start_index}"
+ slice_name = f"i,{start_index},k"
+ else:
+ start_index = slice_range[0] + 1
+ end_index = slice_range[1]
+ slice_title = f", slice j={start_index}:{end_index}"
+ slice_name = f"i,{start_index}:{end_index},k"
+ xc, yc = get_xz_coords(cfg, data, n)
+ mx = 2 * data.nx - 1
+ my = 2 * data.nz - 1
+ xname = "x"
+ yname = "z"
+ return xc, yc, slice_title, slice_name, mx, my, xname, yname
+
+
+
+
+[docs]
+def get_xy_slice(
+ cfg: PlopmConfig, data: SimData, n: int
+) -> tuple[NDArray, NDArray, str, str, int, int, str, str]:
+ """Prepare geometry and labels for an xy slice.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map and slice configuration.
+ data : SimData
+ Loaded grid data.
+ n : int
+ Slice index.
+
+ Returns
+ -------
+ tuple
+ Coordinate meshes, display and filename slice labels, mapped grid
+ dimensions, and coordinate-axis names.
+
+ """
+ slice_range = cfg.slice[n][2]
+ nz = data.nz
+ if slice_range[0] == ":":
+ cfg.slice[n][2] = [0, nz]
+ slice_title = f", slice k={1}:{nz}"
+ slice_name = f"i,j,{1}:{nz}"
+ elif slice_range[0] == slice_range[1] - 1:
+ start_index = slice_range[0] + 1
+ slice_title = f", slice k={start_index}"
+ slice_name = f"i,j,{start_index}"
+ else:
+ start_index = slice_range[0] + 1
+ end_index = slice_range[1]
+ slice_title = f", slice k={start_index}:{end_index}"
+ slice_name = f"i,j,{start_index}:{end_index}"
+ xc, yc = get_xy_coords(cfg, data, n)
+ mx = 2 * data.nx - 1
+ my = 2 * data.ny - 1
+ xname = "x"
+ yname = "y"
+ return xc, yc, slice_title, slice_name, mx, my, xname, yname
+
+
+
+
+[docs]
+def transform_grid(
+ cfg: PlopmConfig, n: int, xc: NDArray, yc: NDArray
+) -> tuple[NDArray, NDArray]:
+ """Rotate and translate a two-dimensional coordinate mesh.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Rotation and translation settings.
+ n : int
+ Map index.
+ xc, yc : np.ndarray
+ Coordinate meshes to transform.
+
+ Returns
+ -------
+ tuple[np.ndarray, np.ndarray]
+ Transformed x- and y-coordinate meshes.
+
+ """
+ grd = int(cfg.rotation[n])
+ angle = grd * np.pi / 180
+ cos_val = np.cos(angle)
+ sin_val = np.sin(angle)
+ length = xc[-1][-1] - xc[0][0]
+ width = yc[0][-1] - yc[-1][0]
+ x_dis = float(cfg.translation[n][0][1:])
+ y_dis = float(cfg.translation[n][1][:-1])
+ base_x = 1.5 * length
+ base_y = 1.5 * width
+ dx = xc - base_x
+ dy = yc - base_y
+ return (
+ base_x + x_dis + dx * cos_val - dy * sin_val,
+ base_y + y_dis + dy * cos_val + dx * sin_val,
+ )
+
+
+
+
+[docs]
+def map_xz(
+ cfg: PlopmConfig,
+ data: SimData,
+ var: str,
+ values: NDArray,
+ n: int,
+ mx: int,
+ my: int,
+ features: list | None = None,
+ feature_id: int = 1,
+) -> NDArray:
+ """Aggregate active-cell values onto an xz slice.
+
+ Values are aggregated through the selected j interval. Permeability uses
+ arithmetic or harmonic thickness weighting according to flow direction;
+ other properties use the configured aggregation or pore-volume weighting.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Slice and aggregation configuration.
+ data : SimData
+ Loaded grid properties and active-cell mapping.
+ var : str
+ Variable name.
+ values : np.ndarray
+ Values in active-cell order.
+ n : int
+ Map index.
+ mx, my : int
+ Mapped grid dimensions.
+ features : list, optional
+ Wells or faults grouped by label.
+ feature_id : int, default: 1
+ Category assigned when mapping one feature.
+
+ Returns
+ -------
+ np.ndarray
+ Values on the flattened xz plotting grid.
+
+ """
+ how = cfg.aggregation[n]
+ nx = data.nx
+ ny = data.ny
+ nz = data.nz
+ slice_start, slice_end = cfg.slice[n][1]
+ layer_size = nx * ny
+ porv = data.porv
+ active_idx = data.active_idx
+ dy = data.dy
+ mapped_values = np.full(mx * my, np.nan)
+ is_wells_or_faults = features is not None
+ is_sum_property = var in cfg.mass_vars or var in [
+ "porv",
+ "dy",
+ "tranx",
+ "tranz",
+ ]
+ is_caprock = var in cfg.caprock_vars
+ is_arithmetic_perm = var in ["permx", "permz"]
+ for k in range(nz):
+ layer_offset = k * layer_size
+ output_layer_offset = 2 * (nz - k - 1) * mx
+ for i in range(nx):
+ p_v, val, d_y = 0.0, 0.0, 0.0
+ if how == "min":
+ val = np.inf
+ if how == "max":
+ val = -np.inf
+ for sld in range(slice_start, slice_end):
+ ind = i + sld * nx + layer_offset
+ cell_pv = porv[ind]
+ if cell_pv > 0:
+ active_id = active_idx[ind]
+ if how and not is_wells_or_faults:
+ if how == "first":
+ p_v = 1.0
+ if var == "index_i":
+ val = i + 1
+ elif var == "index_j":
+ val = sld + 1
+ elif var == "index_k":
+ val = k + 1
+ else:
+ val = values[active_id]
+ break
+ if how == "last":
+ p_v = 1.0
+ if var == "index_i":
+ val = i + 1
+ elif var == "index_j":
+ val = sld + 1
+ elif var == "index_k":
+ val = k + 1
+ else:
+ val = values[active_id]
+ elif how == "min":
+ p_v = 1.0
+ val = min(val, values[active_id])
+ elif how == "max":
+ p_v = 1.0
+ val = max(val, values[active_id])
+ elif how == "sum":
+ p_v = 1.0
+ val += values[active_id]
+ elif how == "mean":
+ p_v += 1.0
+ val += values[active_id]
+ elif how == "pvmean":
+ p_v += cell_pv
+ val += values[active_id] * cell_pv
+ elif how == "harmonic":
+ cell_value = values[active_id]
+ d_y += dy[active_id]
+ val = (
+ np.inf
+ if cell_value == 0
+ else val + dy[active_id] / cell_value
+ )
+ p_v += cell_pv
+ elif how == "arithmetic":
+ p_v += dy[active_id]
+ val += values[active_id] * dy[active_id]
+ elif is_sum_property:
+ p_v = 1.0
+ val += values[active_id]
+ elif is_caprock:
+ p_v = 1.0
+ val = values[active_id]
+ break
+ elif is_arithmetic_perm:
+ p_v += dy[active_id]
+ val += values[active_id] * dy[active_id]
+ elif var == "permy":
+ cell_value = values[active_id]
+ p_v = 1
+ d_y += dy[active_id]
+ val = (
+ np.inf
+ if cell_value == 0
+ else val + dy[active_id] / cell_value
+ )
+ elif var == "grid":
+ p_v = 1
+ val = 1
+ elif var in ["wells", "faults"]:
+ p_v = 1
+ val = feature_id
+ elif var == "index_i":
+ p_v = 1
+ val = i + 1
+ elif var == "index_j":
+ p_v = 1
+ val = sld + 1
+ elif var == "index_k":
+ p_v = 1
+ val = k + 1
+ else:
+ p_v += cell_pv
+ val += values[active_id] * cell_pv
+ if how == "harmonic" or (not how and var == "permy"):
+ mapped_values[2 * i + output_layer_offset] = (
+ np.nan
+ if p_v == 0
+ else 0.0 if val == np.inf else np.nan if val == 0 else d_y / val
+ )
+ else:
+ mapped_values[2 * i + output_layer_offset] = (
+ np.nan if p_v == 0 else val / p_v
+ )
+ if is_wells_or_faults:
+ assert features is not None
+ for index, vals in enumerate(features):
+ for value in vals:
+ if value:
+ for k in range(value[2], value[3] + 1):
+ ind = value[0] + value[1] * nx + k * layer_size
+ if not cfg.global_range:
+ if porv[ind] > 0 and slice_start <= value[1] < slice_end:
+ mapped_values[2 * value[0] + 2 * (nz - k - 1) * mx] = (
+ index + 1
+ )
+ else:
+ if porv[ind] > 0:
+ mapped_values[2 * value[0] + 2 * (nz - k - 1) * mx] = (
+ index + 1
+ )
+ return mapped_values
+
+
+
+
+[docs]
+def map_yz(
+ cfg: PlopmConfig,
+ data: SimData,
+ var: str,
+ values: NDArray,
+ n: int,
+ mx: int,
+ my: int,
+ features: list | None = None,
+ feature_id: int = 1,
+) -> NDArray:
+ """Aggregate active-cell values onto a yz slice.
+
+ Values are aggregated through the selected i interval. Permeability uses
+ arithmetic or harmonic thickness weighting according to flow direction;
+ other properties use the configured aggregation or pore-volume weighting.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Slice and aggregation configuration.
+ data : SimData
+ Loaded grid properties and active-cell mapping.
+ var : str
+ Variable name.
+ values : np.ndarray
+ Values in active-cell order.
+ n : int
+ Map index.
+ mx, my : int
+ Mapped grid dimensions.
+ features : list, optional
+ Wells or faults grouped by label.
+ feature_id : int, default: 1
+ Category assigned when mapping one feature.
+
+ Returns
+ -------
+ np.ndarray
+ Values on the flattened yz plotting grid.
+
+ """
+ how = cfg.aggregation[n]
+ nx = data.nx
+ ny = data.ny
+ nz = data.nz
+ slice_start, slice_end = cfg.slice[n][0]
+ layer_size = nx * ny
+ porv = data.porv
+ active_idx = data.active_idx
+ dx = data.dx
+ mapped_values = np.full(mx * my, np.nan)
+ is_wells_or_faults = features is not None
+ is_sum_property = var in cfg.mass_vars or var in [
+ "porv",
+ "dx",
+ "trany",
+ "tranz",
+ ]
+ is_caprock = var in cfg.caprock_vars
+ is_arithmetic_perm = var in ["permy", "permz"]
+ for k in range(nz):
+ layer_offset = k * layer_size
+ output_layer_offset = 2 * (nz - k - 1) * mx
+ for j in range(ny):
+ row_offset = j * nx
+ p_v, val, d_x = 0.0, 0.0, 0.0
+ if how == "min":
+ val = np.inf
+ if how == "max":
+ val = -np.inf
+ for sld in range(slice_start, slice_end):
+ ind = sld + row_offset + layer_offset
+ cell_pv = porv[ind]
+ if cell_pv > 0:
+ active_id = active_idx[ind]
+ if how and not is_wells_or_faults:
+ if how == "first":
+ p_v = 1.0
+ if var == "index_i":
+ val = sld + 1
+ elif var == "index_j":
+ val = j + 1
+ elif var == "index_k":
+ val = k + 1
+ else:
+ val = values[active_id]
+ break
+ if how == "last":
+ p_v = 1.0
+ if var == "index_i":
+ val = sld + 1
+ elif var == "index_j":
+ val = j + 1
+ elif var == "index_k":
+ val = k + 1
+ else:
+ val = values[active_id]
+ elif how == "min":
+ p_v = 1.0
+ val = min(val, values[active_id])
+ elif how == "max":
+ p_v = 1.0
+ val = max(val, values[active_id])
+ elif how == "sum":
+ p_v = 1.0
+ val += values[active_id]
+ elif how == "mean":
+ p_v += 1.0
+ val += values[active_id]
+ elif how == "pvmean":
+ p_v += cell_pv
+ val += values[active_id] * cell_pv
+ elif how == "harmonic":
+ cell_value = values[active_id]
+ d_x += dx[active_id]
+ val = (
+ np.inf
+ if cell_value == 0
+ else val + dx[active_id] / cell_value
+ )
+ p_v += cell_pv
+ elif how == "arithmetic":
+ p_v += dx[active_id]
+ val += values[active_id] * dx[active_id]
+ elif is_sum_property:
+ p_v = 1.0
+ val += values[active_id]
+ elif is_caprock:
+ p_v = 1.0
+ val = values[active_id]
+ break
+ elif is_arithmetic_perm:
+ p_v += dx[active_id]
+ val += values[active_id] * dx[active_id]
+ elif var == "permx":
+ cell_value = values[active_id]
+ p_v = 1
+ d_x += dx[active_id]
+ val = (
+ np.inf
+ if cell_value == 0
+ else val + dx[active_id] / cell_value
+ )
+ elif var == "grid":
+ p_v = 1
+ val = 1
+ elif var in ["wells", "faults"]:
+ p_v = 1
+ val = feature_id
+ elif var == "index_i":
+ p_v = 1
+ val = sld + 1
+ elif var == "index_j":
+ p_v = 1
+ val = j + 1
+ elif var == "index_k":
+ p_v = 1
+ val = k + 1
+ else:
+ p_v += cell_pv
+ val += values[active_id] * cell_pv
+ if how == "harmonic" or (not how and var == "permx"):
+ mapped_values[2 * j + output_layer_offset] = (
+ np.nan
+ if p_v == 0
+ else 0.0 if val == np.inf else np.nan if val == 0 else d_x / val
+ )
+ else:
+ mapped_values[2 * j + output_layer_offset] = (
+ np.nan if p_v == 0 else val / p_v
+ )
+ if is_wells_or_faults:
+ assert features is not None
+ for index, vals in enumerate(features):
+ for value in vals:
+ if value:
+ for k in range(value[2], value[3] + 1):
+ ind = value[0] + value[1] * nx + k * layer_size
+ if not cfg.global_range:
+ if porv[ind] > 0 and slice_start <= value[0] < slice_end:
+ mapped_values[2 * value[1] + 2 * (nz - k - 1) * mx] = (
+ index + 1
+ )
+ else:
+ if porv[ind] > 0:
+ mapped_values[2 * value[1] + 2 * (nz - k - 1) * mx] = (
+ index + 1
+ )
+ return mapped_values
+
+
+
+
+[docs]
+def map_xy(
+ cfg: PlopmConfig,
+ data: SimData,
+ var: str,
+ values: NDArray,
+ n: int,
+ mx: int,
+ my: int,
+ features: list | None = None,
+ feature_id: int = 1,
+) -> NDArray:
+ """Aggregate active-cell values onto an xy slice.
+
+ Values are aggregated through the selected k interval. Dual-porosity rows
+ are included when enabled, and permeability is weighted according to the
+ vertical flow direction.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Slice, aggregation, and dual-grid configuration.
+ data : SimData
+ Loaded grid properties and active-cell mapping.
+ var : str
+ Variable name.
+ values : np.ndarray
+ Values in active-cell order.
+ n : int
+ Map index.
+ mx, my : int
+ Mapped grid dimensions.
+ features : list, optional
+ Wells or faults grouped by label.
+ feature_id : int, default: 1
+ Category assigned when mapping one feature.
+
+ Returns
+ -------
+ np.ndarray
+ Values on the flattened xy plotting grid.
+
+ """
+ how = cfg.aggregation[n]
+ nx = data.nx
+ ny_total = data.ny
+ dual = cfg.dual_grid[n] == "1" if n < len(cfg.dual_grid) else False
+ ny = int((ny_total - 1) / 2) if dual else ny_total
+ slice_start, slice_end = cfg.slice[n][2]
+ layer_size = nx * ny_total
+ porv = data.porv
+ active_idx = data.active_idx
+ dz = data.dz
+ mapped_values = np.full(mx * my, np.nan)
+ is_wells_or_faults = features is not None
+ is_sum_property = var in cfg.mass_vars or var in [
+ "porv",
+ "dz",
+ "tranx",
+ "trany",
+ ]
+ is_caprock = var in cfg.caprock_vars
+ is_arithmetic_perm = var in ["permx", "permy"]
+ for j in range(ny):
+ row_offset = j * nx
+ dual_row_offset = (j + ny + 1) * nx
+ for i in range(nx):
+ p_v, val, d_z = 0.0, 0.0, 0.0
+ if how == "min":
+ val = np.inf
+ if how == "max":
+ val = -np.inf
+ for sld in range(slice_start, slice_end):
+ layer_offset = sld * layer_size
+ ind = i + row_offset + layer_offset
+ idd = i + dual_row_offset + layer_offset
+ cell_pv = porv[ind]
+ dual_cell_pv = porv[idd] if dual else 0
+ if cell_pv > 0 or (dual and dual_cell_pv > 0):
+ active_id = active_idx[ind]
+ dual_active_id = active_idx[idd] if dual else active_id
+ if how and not is_wells_or_faults:
+ if how == "first":
+ p_v = 1.0
+ if var == "index_i":
+ val = i + 1
+ elif var == "index_j":
+ val = j + 1
+ elif var == "index_k":
+ val = sld + 1
+ else:
+ val = values[active_id]
+ break
+ if how == "last":
+ p_v = 1.0
+ if var == "index_i":
+ val = i + 1
+ elif var == "index_j":
+ val = j + 1
+ elif var == "index_k":
+ val = sld + 1
+ else:
+ val = values[active_id]
+ elif how == "min":
+ p_v = 1.0
+ if cell_pv > 0:
+ val = min(val, values[active_id])
+ if dual and dual_cell_pv > 0:
+ val = min(val, values[dual_active_id])
+ elif how == "max":
+ p_v = 1.0
+ if cell_pv > 0:
+ val = max(val, values[active_id])
+ if dual and dual_cell_pv > 0:
+ val = max(val, values[dual_active_id])
+ elif how == "sum":
+ p_v = 1.0
+ if cell_pv > 0:
+ val += values[active_id]
+ if dual and dual_cell_pv > 0:
+ val += values[dual_active_id]
+ elif how == "mean":
+ if cell_pv > 0:
+ p_v += 1.0
+ val += values[active_id]
+ if dual and dual_cell_pv > 0:
+ p_v += 1.0
+ val += values[dual_active_id]
+ elif how == "pvmean":
+ if cell_pv > 0:
+ p_v += cell_pv
+ val += values[active_id] * cell_pv
+ if dual and dual_cell_pv > 0:
+ p_v += dual_cell_pv
+ val += values[dual_active_id] * dual_cell_pv
+ elif how == "harmonic":
+ if cell_pv > 0:
+ cell_value = values[active_id]
+ d_z += dz[active_id]
+ val = (
+ np.inf
+ if cell_value == 0
+ else val + dz[active_id] / cell_value
+ )
+ p_v += cell_pv
+ if dual and dual_cell_pv > 0:
+ cell_value = values[dual_active_id]
+ d_z += dz[dual_active_id]
+ val = (
+ np.inf
+ if cell_value == 0
+ else val + dz[dual_active_id] / cell_value
+ )
+ p_v += dual_cell_pv
+ elif how == "arithmetic":
+ if cell_pv > 0:
+ p_v += dz[active_id]
+ val += values[active_id] * dz[active_id]
+ if dual and dual_cell_pv > 0:
+ p_v += dz[dual_active_id]
+ val += values[dual_active_id] * dz[dual_active_id]
+ elif is_sum_property:
+ p_v = 1.0
+ if cell_pv > 0:
+ val += values[active_id]
+ if dual and dual_cell_pv > 0:
+ val += values[dual_active_id]
+ elif is_caprock:
+ p_v = 1.0
+ val = values[active_id]
+ break
+ elif is_arithmetic_perm:
+ if cell_pv > 0:
+ p_v += dz[active_id]
+ val += values[active_id] * dz[active_id]
+ if dual and dual_cell_pv > 0:
+ p_v += dz[dual_active_id]
+ val += values[dual_active_id] * dz[dual_active_id]
+ elif var == "permz":
+ p_v = 1
+ if cell_pv > 0:
+ cell_value = values[active_id]
+ d_z += dz[active_id]
+ val = (
+ np.inf
+ if cell_value == 0
+ else val + dz[active_id] / cell_value
+ )
+ if dual and dual_cell_pv > 0:
+ cell_value = values[dual_active_id]
+ d_z += dz[dual_active_id]
+ val = (
+ np.inf
+ if cell_value == 0
+ else val + dz[dual_active_id] / cell_value
+ )
+ elif var == "grid":
+ p_v = 1
+ val = 1
+ elif var in ["wells", "faults"]:
+ p_v = 1
+ val = feature_id
+ elif var == "index_i":
+ p_v = 1
+ val = i + 1
+ elif var == "index_j":
+ p_v = 1
+ val = j + 1
+ elif var == "index_k":
+ p_v = 1
+ val = sld + 1
+ else:
+ if cell_pv > 0:
+ p_v += cell_pv
+ val += values[active_id] * cell_pv
+ if dual and dual_cell_pv > 0:
+ p_v += dual_cell_pv
+ val += values[dual_active_id] * dual_cell_pv
+ if how == "harmonic" or (not how and var == "permz"):
+ mapped_values[2 * i + 2 * j * mx] = (
+ np.nan
+ if p_v == 0
+ else 0.0 if val == np.inf else np.nan if val == 0 else d_z / val
+ )
+ else:
+ mapped_values[2 * i + 2 * j * mx] = np.nan if p_v == 0 else val / p_v
+ if is_wells_or_faults:
+ assert features is not None
+ for index, vals in enumerate(features):
+ for value in vals:
+ if value:
+ for k in range(value[2], value[3] + 1):
+ ind = value[0] + value[1] * nx + k * layer_size
+ if not cfg.global_range:
+ if porv[ind] > 0 and slice_start <= k < slice_end:
+ mapped_values[2 * value[0] + 2 * value[1] * mx] = (
+ index + 1
+ )
+ else:
+ if porv[ind] > 0:
+ mapped_values[2 * value[0] + 2 * value[1] * mx] = (
+ index + 1
+ )
+ if dual and cfg.difference_input:
+ mapped_values = mapped_values[: (2 * nx - 1) * (2 * ny - 1)]
+ return mapped_values
+
+
+
+
+
+
+# SPDX-FileCopyrightText: 2024-2026 NORCE Research AS
+# SPDX-License-Identifier: GPL-3.0
+# pylint: disable=R0911,R0912,R0913,R0915,R0917,R1702,R0914,C0302,E1102
+
+"""Read and derive plotting quantities from OPM Flow output.
+
+The module opens INIT, UNRST, EGRID, SMSPEC, deck, and CSV data; constructs
+plotting coordinates; evaluates variable expressions; and derives saturation,
+mass, caprock, distance, well, and fault quantities.
+"""
+
+import csv
+import datetime
+import os
+import sys
+from contextlib import nullcontext
+
+import numpy as np
+from alive_progress import alive_bar
+from numpy.typing import NDArray
+from opm.io.ecl import EclFile as OpmFile
+from opm.io.ecl import EGrid as OpmGrid
+from opm.io.ecl import ERst as OpmRestart
+from opm.io.ecl import ESmry as OpmSummary
+
+from plopm.config.config import PlopmConfig, SimData
+from plopm.utils.initialization import mass_unit, spatial_unit
+from plopm.utils.terminal import (
+ cli_error_value,
+ cli_info_value,
+ plopm_error,
+ plopm_info,
+)
+
+csv.field_size_limit(sys.maxsize)
+
+GAS_DEN_REF = 1.86843
+WAT_DEN_REF = 998.108
+
+
+
+[docs]
+def read_case(
+ deck: str,
+ gif: bool,
+ vtk: bool,
+ variables: list,
+ restart: list,
+ filters: list,
+ n: int = 0,
+) -> SimData:
+ """Open the OPM output required for one simulation case.
+
+ Parameters
+ ----------
+ deck : str
+ Simulation-case stem without an extension.
+ gif, vtk : bool
+ Output modes controlling restart and grid loading.
+ variables : list
+ Requested variables or expressions.
+ restart : list
+ Requested restart report steps.
+ filters : list
+ Property-filter expressions.
+ n : int, default: 0
+ Case index used to select per-case settings.
+
+ Returns
+ -------
+ SimData
+ Loaded readers, grid properties, and report-step metadata.
+
+ """
+ if os.path.isfile(f"{deck}.INIT"):
+ init = OpmFile(f"{deck}.INIT")
+ else:
+ plopm_error(f"unable to find {cli_error_value(f'{deck}.INIT')}")
+ unrst = OpmRestart(f"{deck}.UNRST") if os.path.isfile(f"{deck}.UNRST") else None
+ egrid = (
+ OpmGrid(f"{deck}.EGRID")
+ if os.path.isfile(f"{deck}.EGRID") and not vtk
+ else None
+ )
+
+ porv = np.array(init["PORV"])
+ dx = np.array(init["DX"])
+ dy = np.array(init["DY"])
+ dz = np.array(init["DZ"])
+
+ act_mask = porv > 0
+ pv = porv[act_mask]
+ actind = np.cumsum(act_mask) - 1
+
+ tnrst = []
+ ntot = 1
+
+ if filters[n]:
+ porv0 = porv.copy()
+ for value in filters[n].split("&"):
+ filte = value.strip().split(" ")
+ key = filte[0].upper()
+ if init.count(key):
+ arr = np.array(init[key])
+ mask = porv0 > 0
+ porv[mask] = _apply_filter(porv[mask], arr, filte[1], float(filte[2]))
+
+ if unrst:
+ steps = unrst.report_steps
+ ntot = steps[-1] + 1
+ if unrst.count("DOUBHEAD", 0):
+ tnrst = [unrst["DOUBHEAD", ntm][0] for ntm in steps]
+ if restart[0] == -1:
+ restart = unrst.report_steps if gif else [ntot - 1]
+ elif restart[0] == -1:
+ restart = [ntot - 1]
+
+ nx = ny = nz = 0
+
+ if egrid:
+ dim = egrid.dimension
+ nx, ny, nz = dim
+ elif "index_i" in variables or "index_j" in variables or "index_k" in variables:
+ grid = OpmGrid(f"{deck}.EGRID")
+ dim = grid.dimension
+ nx, ny, nz = dim
+
+ if not tnrst:
+ tnrst = [0] * len(restart)
+
+ return SimData(
+ init,
+ unrst,
+ egrid,
+ porv,
+ dx,
+ dy,
+ dz,
+ pv,
+ actind,
+ restart,
+ tnrst,
+ porv.size,
+ ntot,
+ nx,
+ ny,
+ nz,
+ )
+
+
+
+
+[docs]
+def get_yz_coords(cfg: PlopmConfig, data: SimData, n: int) -> tuple[NDArray, NDArray]:
+ """Build coordinate meshes for a yz slice.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map configuration.
+ data : SimData
+ Loaded grid data.
+ n : int
+ Slice index.
+
+ Returns
+ -------
+ tuple[np.ndarray, np.ndarray]
+ Y- and z-coordinate meshes.
+
+ """
+ xyz_func = data.grid.xyz_from_ijk
+ ny_val = data.ny
+ nz_val = data.nz
+ base_i_all = cfg.slice[n][0][0]
+ total_size = nz_val * 4 * ny_val
+ xc_list = [0] * total_size
+ yc_list = [0] * total_size
+ idx = 0
+ for j in range(nz_val):
+ base_k = nz_val - j - 1
+ base_idx_second = idx + 2 * ny_val
+ tmp_idx = base_idx_second
+ for i in range(ny_val):
+ val = xyz_func(base_i_all, i, base_k, True)
+ xc_list[idx] = val[1][4]
+ yc_list[idx] = val[2][4]
+ idx += 1
+ xc_list[idx] = val[1][6]
+ yc_list[idx] = val[2][6]
+ idx += 1
+ xc_list[tmp_idx] = val[1][0]
+ yc_list[tmp_idx] = val[2][0]
+ tmp_idx += 1
+ xc_list[tmp_idx] = val[1][2]
+ yc_list[tmp_idx] = val[2][2]
+ tmp_idx += 1
+ idx = base_idx_second + 2 * ny_val
+ xc_array = np.asarray(xc_list)
+ yc_array = np.asarray(yc_list)
+ return xc_array.reshape(2 * nz_val, 2 * ny_val), yc_array.reshape(
+ 2 * nz_val, 2 * ny_val
+ )
+
+
+
+
+[docs]
+def get_xz_coords(cfg: PlopmConfig, data: SimData, n: int) -> tuple[NDArray, NDArray]:
+ """Build coordinate meshes for an xz slice.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map configuration.
+ data : SimData
+ Loaded grid data.
+ n : int
+ Slice index.
+
+ Returns
+ -------
+ tuple[np.ndarray, np.ndarray]
+ X- and z-coordinate meshes.
+
+ """
+ xyz_func = data.grid.xyz_from_ijk
+ nx_val = data.nx
+ nz_val = data.nz
+ base_j_all = cfg.slice[n][1][0]
+ total_size = nz_val * 4 * nx_val
+ xc_list = [0] * total_size
+ yc_list = [0] * total_size
+ idx = 0
+ for j in range(nz_val):
+ base_k = nz_val - j - 1
+ base_idx_second = idx + 2 * nx_val
+ tmp_idx = base_idx_second
+ for i in range(nx_val):
+ val = xyz_func(i, base_j_all, base_k, True)
+ xc_list[idx] = val[0][4]
+ yc_list[idx] = val[2][4]
+ idx += 1
+ xc_list[idx] = val[0][5]
+ yc_list[idx] = val[2][5]
+ idx += 1
+ xc_list[tmp_idx] = val[0][0]
+ yc_list[tmp_idx] = val[2][0]
+ tmp_idx += 1
+ xc_list[tmp_idx] = val[0][1]
+ yc_list[tmp_idx] = val[2][1]
+ tmp_idx += 1
+ idx = base_idx_second + 2 * nx_val
+ xc_array = np.asarray(xc_list)
+ yc_array = np.asarray(yc_list)
+ return xc_array.reshape(2 * nz_val, 2 * nx_val), yc_array.reshape(
+ 2 * nz_val, 2 * nx_val
+ )
+
+
+
+
+[docs]
+def get_xy_coords(cfg: PlopmConfig, data: SimData, n: int) -> tuple[NDArray, NDArray]:
+ """Build coordinate meshes for an xy slice.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map configuration.
+ data : SimData
+ Loaded grid data.
+ n : int
+ Slice index.
+
+ Returns
+ -------
+ tuple[np.ndarray, np.ndarray]
+ X- and y-coordinate meshes.
+
+ """
+ xyz_func = data.grid.xyz_from_ijk
+ nx_val = data.nx
+ ny_val = data.ny
+ base_k_all = cfg.slice[n][2][0]
+ total_size = ny_val * 4 * nx_val
+ xc_list = [0] * total_size
+ yc_list = [0] * total_size
+ idx = 0
+ for j in range(ny_val):
+ base_idx_second = idx + 2 * nx_val
+ tmp_idx = base_idx_second
+ for i in range(nx_val):
+ val = xyz_func(i, j, base_k_all, True)
+ xc_list[idx] = val[0][0]
+ yc_list[idx] = val[1][0]
+ idx += 1
+ xc_list[idx] = val[0][1]
+ yc_list[idx] = val[1][1]
+ idx += 1
+ xc_list[tmp_idx] = val[0][2]
+ yc_list[tmp_idx] = val[1][2]
+ tmp_idx += 1
+ xc_list[tmp_idx] = val[0][3]
+ yc_list[tmp_idx] = val[1][3]
+ tmp_idx += 1
+ idx = base_idx_second + 2 * nx_val
+ xc_array = np.asarray(xc_list)
+ yc_array = np.asarray(yc_list)
+ return xc_array.reshape(2 * ny_val, 2 * nx_val), yc_array.reshape(
+ 2 * ny_val, 2 * nx_val
+ )
+
+
+
+def _resolve_var(
+ cfg: PlopmConfig,
+ data: SimData,
+ key_up: str,
+ key_low: str,
+ step: int,
+ init: OpmFile,
+ unrst: OpmRestart,
+ mass_all: list,
+ caprock_list: list,
+):
+ """Resolve a variable from stored or derived quantities.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized configuration.
+ data : SimData
+ Loaded simulation data.
+ key_up, key_low : str
+ OPM keyword and normalized variable name.
+ step : int
+ Restart report step.
+ init : OpmFile
+ INIT reader.
+ unrst : OpmRestart
+ UNRST reader.
+ mass_all, caprock_list : list
+ Supported derived variable names.
+
+ Returns
+ -------
+ np.ndarray or None
+ Resolved values, or ``None`` when unavailable.
+
+ """
+ if init.count(key_up):
+ return 1.0 * init[key_up, 0]
+ if unrst is not None and unrst.count(key_up, step):
+ return 1.0 * unrst[key_up, step]
+ if key_low in mass_all:
+ return _get_mass(data, key_low, step)
+ if key_low in caprock_list:
+ val, _ = _get_caprock(data, key_low, step, cfg.stress_coefficient)
+ return val
+ if key_low in ["swat", "soil", "sgas"]:
+ return _get_saturation(data.unrst, key_low, step)
+ return None
+
+
+def _read_histogram(
+ cfg: PlopmConfig, data: SimData, tokens: list, step: int
+) -> NDArray:
+ """Read values used to create a histogram.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized configuration.
+ data : SimData
+ Loaded simulation data.
+ tokens : list
+ Parsed variable-expression tokens.
+ step : int
+ Restart report step.
+
+ Returns
+ -------
+ np.ndarray
+ Values in global cell order with inactive cells set to NaN.
+
+ """
+ quan0_low = tokens[0]
+ quan0 = quan0_low.upper()
+ porv = data.porv
+ nxyz = data.ncells
+ init = data.init
+ unrst = data.unrst
+ mass_all = cfg.mass_vars + cfg.mass_fracs
+ caprock_list = cfg.caprock_vars
+ if quan0 != "PORV":
+ act = porv > 0
+ else:
+ act = porv > -1
+ var = np.nan * np.ones(nxyz, dtype=float)
+ result = _resolve_var(
+ cfg, data, quan0, quan0_low, step, init, unrst, mass_all, caprock_list
+ )
+ if result is not None:
+ var[act] = result
+ else:
+ plopm_error(f"not found {cli_error_value(f'-v {tokens[0]}')}.")
+ if len(tokens) > 1:
+ ops = tokens[1::2]
+ for j, val in enumerate(tokens[2::2]):
+ val_up = val.upper()
+ if val[0].isdigit() and not val[-1].isdigit():
+ if unrst is None:
+ plopm_error(f"not found {cli_error_value(f'-v {val}')}.")
+ other = 1.0 * unrst[val[1:].upper(), int(val[0])]
+ elif val[0].isdigit() and val[-1].isdigit():
+ other = np.full_like(var[act], float(val))
+ else:
+ other = _resolve_var(
+ cfg,
+ data,
+ val_up,
+ val,
+ step,
+ init,
+ unrst,
+ mass_all,
+ caprock_list,
+ )
+ if other is None:
+ plopm_error(f"not found {cli_error_value(f'-v {val}')}.")
+ var_act = var[act]
+ var[act] = _apply_operator(var_act, other, ops[j])
+ return var
+
+
+def _compute_distance(
+ cfg: PlopmConfig, data: SimData, tokens: list, n: int
+) -> tuple[NDArray, NDArray]:
+ """Compute distance from selected cells to target points.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Distance and sensor configuration.
+ data : SimData
+ Loaded simulation data.
+ tokens : list
+ Parsed expression selecting active cells.
+ n : int
+ Plot index.
+
+ Returns
+ -------
+ tuple[np.ndarray, np.ndarray]
+ Finite distances and their simulation times.
+
+ """
+ xyz_func = data.grid.xyz_from_ijk
+ nx_val = data.nx
+ ny_val = data.ny
+ nz_val = data.nz
+ nxyz = data.ncells
+ ntot = data.nsteps
+ porv = data.porv
+ init = data.init
+ unrst = data.unrst
+ mass_all = cfg.mass_vars + cfg.mass_fracs
+ distance_type = cfg.distance[0]
+ xyz = np.zeros((nxyz, 3), dtype=float)
+ act = porv > 0
+ time = np.array(data.times)
+ distance = np.nan * np.ones(ntot)
+ index = 0
+ for k in range(nz_val):
+ for j in range(ny_val):
+ for i in range(nx_val):
+ xyz[index, :] = np.mean(xyz_func(i, j, k, True), axis=1)
+ index += 1
+ if cfg.distance[1] == "sensor":
+ ind = (
+ cfg.slice[n][0]
+ + cfg.slice[n][1] * nx_val
+ + cfg.slice[n][2] * nx_val * ny_val
+ )
+ points = [xyz[ind, :]]
+ sensor_loc = f"[{points[0][0]:.2E},{points[0][1]:.2E},{points[0][2]:.2E}]"
+ plopm_info(
+ f"computing the {cli_info_value(cfg.distance[0])} distance of "
+ f"{cli_info_value(tokens[0])} to the sensor "
+ f"{cli_info_value(sensor_loc)} [m]"
+ )
+ else:
+ points = []
+ for k in range(nz_val):
+ if ny_val > 1:
+ base_k = k * nx_val * ny_val
+ for i in range(nx_val):
+ ind = i + base_k
+ if act[ind]:
+ points.append(xyz[ind])
+ ind = i + (ny_val - 1) * nx_val + base_k
+ if act[ind]:
+ points.append(xyz[ind])
+ if nx_val > 1:
+ base_k = k * nx_val * ny_val
+ for j in range(ny_val):
+ ind = j * nx_val + base_k
+ if act[ind]:
+ points.append(xyz[ind])
+ ind = nx_val - 1 + j * nx_val + base_k
+ if act[ind]:
+ points.append(xyz[ind])
+ plopm_info(
+ f"computing the {cli_info_value(cfg.distance[0])} distance of "
+ f"{cli_info_value(tokens[0])} to the model boundaries"
+ )
+ show_progress = sys.stdout.isatty()
+ if show_progress:
+ bar_ctx = alive_bar(ntot * len(points), bar="fish")
+ else:
+ bar_ctx = nullcontext()
+ with bar_ctx as bar_animation:
+ for step in unrst.report_steps:
+ xyzt = np.copy(xyz)
+ var = np.nan * np.ones(nxyz, dtype=float)
+ quan0_low = tokens[0]
+ quan0_up = tokens[0].upper()
+ if quan0_low in ["index_i", "index_j", "index_k"]:
+ var[act] = _grid_indices(quan0_low, nx_val, ny_val, nz_val)
+ elif unrst.count(quan0_up, step):
+ var[act] = 1.0 * unrst[quan0_up, step]
+ elif quan0_low in mass_all:
+ var[act] = _get_mass(data, quan0_low, step)
+ elif quan0_low in ["swat", "soil", "sgas"]:
+ var[act] = _get_saturation(data.unrst, quan0_low, step)
+ else:
+ flag = f"-dist {','.join(cfg.distance)}"
+ plopm_error(
+ f"invalid {cli_error_value(f'-v {tokens[0]}')} for "
+ f"{cli_info_value(flag)}."
+ )
+ if len(tokens) > 1:
+ ops = tokens[1::2]
+ for j, val in enumerate(tokens[2::2]):
+ val_up = val.upper()
+ if val[0].isdigit() and not val[-1].isdigit():
+ other = 1.0 * unrst[val[1:].upper(), int(val[0])]
+ elif val[0].isdigit() and val[-1].isdigit():
+ other = np.full_like(var[act], float(val))
+ elif init.count(val_up):
+ other = 1.0 * init[val_up, 0]
+ if val_up == "PORV":
+ other = other[act]
+ elif val in ["index_i", "index_j", "index_k"]:
+ var[act] = _grid_indices(val, nx_val, ny_val, nz_val)
+ continue
+ elif unrst.count(val_up, step):
+ other = 1.0 * unrst[val_up, step]
+ elif val in mass_all:
+ other = _get_mass(data, val, step)
+ elif val in ["swat", "soil", "sgas"]:
+ other = _get_saturation(data.unrst, val, step)
+ else:
+ plopm_error(f"not found {cli_error_value(f'-v {val}')}.")
+ var_act = var[act]
+ var[act] = _apply_operator(var_act, other, ops[j])
+ else:
+ var[var > 0] = 1
+ xyzt[var != 1] = np.nan
+ temp = np.nan * np.ones(len(points))
+ for point_index, point in enumerate(points):
+ if show_progress:
+ bar_animation()
+ vals = np.linalg.norm(xyzt - point, axis=1)
+ if not np.all(np.isnan(vals)):
+ if distance_type == "min":
+ temp[point_index] = np.nanmin(vals)
+ else:
+ temp[point_index] = np.nanmax(vals)
+ if not np.isnan(temp).all():
+ if distance_type == "min":
+ distance[step] = np.nanmin(temp)
+ else:
+ distance[step] = np.nanmax(temp)
+ return distance[~np.isnan(distance)], time[~np.isnan(distance)]
+
+
+def _grid_indices(name: str, nx: int, ny: int, nz: int) -> list:
+ """Create one-based grid indices in global cell order.
+
+ Parameters
+ ----------
+ name : {"index_i", "index_j", "index_k"}
+ Grid axis to index.
+ nx, ny, nz : int
+ Grid dimensions.
+
+ Returns
+ -------
+ list[int]
+ One-based indices for all grid cells.
+
+ """
+ nxyz = nx * ny * nz
+ if name == "index_i":
+ return [(grid_index % nx) + 1 for grid_index in range(nxyz)]
+ if name == "index_j":
+ return [((grid_index // nx) % ny) + 1 for grid_index in range(nxyz)]
+ return [(grid_index // (nx * ny)) + 1 for grid_index in range(nxyz)]
+
+
+def _aggregate(var: NDArray, op: str, porv: NDArray) -> NDArray:
+ """_aggregate values with the selected method.
+
+ Parameters
+ ----------
+ var : np.ndarray
+ Values to _aggregate.
+ op : str
+ Aggregation method.
+ porv : np.ndarray
+ Pore-volume weights used by ``"pvmean"``.
+
+ Returns
+ -------
+ np.ndarray or float
+ _aggregated values.
+
+ """
+ if op == "min":
+ return np.min(var)
+ if op == "max":
+ return np.max(var)
+ if op == "sum":
+ return np.sum(var)
+ if op == "mean":
+ return np.mean(var)
+ if op == "pvmean":
+ return np.sum(var * porv) / np.sum(porv)
+ plopm_error(f"unknow/unsupported aggregation {cli_error_value(f'-agg {op}')}.")
+
+
+def _read_values(
+ cfg: PlopmConfig, data: SimData, tokens: list, n: int, ntot: list
+) -> tuple[NDArray, NDArray]:
+ """Read an _aggregated time series or grid-axis profile.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized series configuration.
+ data : SimData
+ Loaded simulation data.
+ tokens : list
+ Parsed variable-expression tokens.
+ n : int
+ Plot index.
+ ntot : list
+ Restart report steps to evaluate.
+
+ Returns
+ -------
+ tuple[np.ndarray, np.ndarray]
+ Values and corresponding time or grid coordinates.
+
+ """
+ slc_val = cfg.slice[n]
+ axis_index = slc_val.index(-1) if -1 in slc_val else -1
+ nx_val = data.nx
+ ny_val = data.ny
+ nz_val = data.nz
+ if axis_index == 0:
+ xsize = nx_val
+ elif axis_index == 1:
+ xsize = ny_val
+ elif axis_index == 2:
+ xsize = nz_val
+ else:
+ xsize = 1
+ if len(ntot) > 1:
+ tsize = len(ntot)
+ time = np.array(data.times)
+ var = 0.0 * np.ones(tsize)
+ else:
+ time = np.array(range(xsize), dtype=float)
+ var = 0.0 * np.ones(xsize)
+ init = data.init
+ unrst = data.unrst
+ mass_all = cfg.mass_vars + cfg.mass_fracs
+ caprock_list = cfg.caprock_vars
+ pv_all = data.active_pv
+ layer_flag = cfg.layer
+ egrid = data.grid
+ quan0_low = tokens[0]
+ quan0_up = quan0_low.upper()
+ ops = tokens[1::2] if len(tokens) > 1 else []
+ for output_index, step in enumerate(ntot):
+ temp = np.ones(xsize, dtype=float)
+ porv = np.ones(xsize, dtype=float)
+ inds = [0] * xsize
+ if layer_flag:
+ if axis_index == 0:
+ for index in range(xsize):
+ inds[index] = egrid.active_index(index, slc_val[1], slc_val[2])
+ elif axis_index == 1:
+ for index in range(xsize):
+ inds[index] = egrid.active_index(slc_val[0], index, slc_val[2])
+ elif axis_index == 2:
+ for index in range(xsize):
+ inds[index] = egrid.active_index(slc_val[0], slc_val[1], index)
+ else:
+ ind0 = egrid.active_index(slc_val[0], slc_val[1], slc_val[2])
+ for index in range(xsize):
+ inds[index] = ind0
+ if quan0_low in mass_all:
+ arr_main = _get_mass(data, quan0_low, step)
+ elif quan0_low in caprock_list:
+ arr_main, _ = _get_caprock(data, quan0_low, step, cfg.stress_coefficient)
+ elif quan0_low in ["swat", "soil", "sgas"]:
+ arr_main = _get_saturation(data.unrst, quan0_low, step)
+ else:
+ arr_main = None
+ if len(tokens) > 1:
+ arr_vals = []
+ for val in tokens[2::2]:
+ if val in mass_all:
+ arr_vals.append(_get_mass(data, val, step))
+ elif val in caprock_list:
+ arr, _ = _get_caprock(data, val, step, cfg.stress_coefficient)
+ arr_vals.append(arr)
+ elif val in ["swat", "soil", "sgas"]:
+ arr_vals.append(_get_saturation(data.unrst, val, step))
+ else:
+ arr_vals.append(np.full_like(temp, np.nan))
+ inds_arr = np.array(inds)
+
+ if unrst.count("RPORV", step):
+ porv = unrst["RPORV", step][inds_arr]
+ else:
+ porv = pv_all[inds_arr]
+
+ if unrst.count(quan0_up, step):
+ temp = 1.0 * unrst[quan0_up, step][inds_arr]
+ # porv-weighted pressure for the dual model
+ if cfg.dual_grid[n] == "1" and cfg.sensor:
+ indd = egrid.active_index(
+ slc_val[0], slc_val[1] + int((data.ny - 1) / 2) + 1, slc_val[2]
+ )
+ presd = unrst[quan0_up, step][indd]
+ if unrst.count("RPORV", step):
+ porvd = unrst["RPORV", step][indd]
+ else:
+ porvd = pv_all[indd]
+ temp = (temp * porv + presd * porvd) / (porv + porvd)
+ elif init.count(quan0_up):
+ temp = 1.0 * init[quan0_up, 0][inds_arr]
+ elif arr_main is not None:
+ temp = arr_main[inds_arr]
+ else:
+ plopm_error(f"not found {cli_error_value(f'-v {tokens[0]}')}.")
+
+ if len(tokens) > 1:
+ for j, val in enumerate(tokens[2::2]):
+ val_up = val.upper()
+ arr_val = arr_vals[j]
+ if val[0].isdigit() and not val[-1].isdigit():
+ other = 1.0 * unrst[val[1:].upper(), int(val[0])][inds_arr]
+ elif val[0].isdigit() and val[-1].isdigit():
+ other = np.full_like(temp, float(val))
+ elif init.count(val_up):
+ other = 1.0 * init[val_up, 0][inds_arr]
+ elif unrst.count(val_up, step):
+ other = 1.0 * unrst[val_up, step][inds_arr]
+ elif not np.isnan(arr_val).all():
+ other = arr_val[inds_arr]
+ else:
+ plopm_error(f"not found {cli_error_value(f'-v {val}')}.")
+ temp = _apply_operator(temp, other, ops[j])
+ ll = np.arange(xsize) + output_index
+ if cfg.aggregation[0]:
+ var[output_index] = _aggregate(temp, cfg.aggregation[0], porv)
+ elif layer_flag:
+ var = temp
+ else:
+ if xsize == 1:
+ var[ll] = temp[0]
+ else:
+ var[ll] = temp
+ if layer_flag and not cfg.aggregation[0]:
+ xyz_func = egrid.xyz_from_ijk
+ if axis_index == 0:
+ for i in range(nx_val):
+ time[i] = np.mean(xyz_func(i, slc_val[1], slc_val[2], True), axis=1)[0]
+ elif axis_index == 1:
+ for j in range(ny_val):
+ time[j] = np.mean(xyz_func(slc_val[0], j, slc_val[2], True), axis=1)[1]
+ else:
+ for k in range(nz_val):
+ time[k] = np.mean(xyz_func(slc_val[0], slc_val[1], k, True), axis=1)[2]
+ return var, time
+
+
+
+[docs]
+def read_series(
+ cfg: PlopmConfig, case: str, values: str, tunit: str, qskl: float, n: int
+) -> tuple[NDArray, NDArray, str, str]:
+ """Read one one-dimensional series for plotting.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized series configuration.
+ case : str
+ Simulation-case stem or CSV path.
+ values : str
+ Variable name or expression.
+ tunit : str
+ Requested time-unit code.
+ qskl : float
+ Scale factor applied to values.
+ n : int
+ Plot or case index.
+
+ Returns
+ -------
+ tuple[np.ndarray, np.ndarray, str, str]
+ Coordinates, values, coordinate unit, and value unit.
+
+ """
+ time, vunit = np.array([0, 1]), ""
+ tskl, tunit = time_unit(tunit)
+ tokens = values.split(" ")
+ csv_flag = cfg.csv_columns[n][0]
+ q0_low = tokens[0]
+ use_sw = "krw" in "".join(cfg.variables)
+ if csv_flag:
+ csvv = np.genfromtxt(f"{case}.csv", delimiter=",", skip_header=1)
+ col_t = cfg.csv_columns[n][0] - 1
+ col_v = cfg.csv_columns[n][1] - 1
+ time = tskl * csvv[:, col_t] / 86400.0
+ var = csvv[:, col_v]
+ elif cfg.distance[0]:
+ xskl, xunit = spatial_unit(cfg.xunits)
+ data = read_case(
+ case, cfg.gif, cfg.vtk, cfg.variables, cfg.restart, cfg.filters
+ )
+ var, time = _compute_distance(cfg, data, tokens, n)
+ vunit = f" ({cfg.distance[0]} distance to {cfg.distance[1]} in {xunit})"
+ var *= xskl
+ elif cfg.histogram[0]:
+ data = read_case(
+ case, cfg.gif, cfg.vtk, cfg.variables, cfg.restart, cfg.filters
+ )
+ var = _read_histogram(cfg, data, tokens, data.steps[0])
+ tunit = ""
+ elif cfg.sensor or cfg.aggregation[0]:
+ data = read_case(
+ case, cfg.gif, cfg.vtk, cfg.variables, cfg.restart, cfg.filters
+ )
+ var, time = _read_values(cfg, data, tokens, n, data.unrst.report_steps)
+ time *= tskl
+ if tunit == "Dates":
+ dates = []
+ unrst = data.unrst
+
+ for step in range(len(unrst)):
+ intehead = unrst["INTEHEAD", step]
+ dates.append(
+ datetime.date(
+ year=int(intehead[66]),
+ month=int(intehead[65]),
+ day=int(intehead[64]),
+ )
+ )
+
+ time = np.asarray(dates)
+ elif cfg.layer:
+ xskl, tunit = spatial_unit(cfg.xunits)
+ data = read_case(
+ case, cfg.gif, cfg.vtk, cfg.variables, cfg.restart, cfg.filters
+ )
+ tmp = data.steps[n] if n < len(cfg.restart) else data.steps[0]
+ var, time = _read_values(cfg, data, tokens, n, [tmp])
+ time *= xskl
+ elif q0_low[:3] in ["krw", "krg"] or q0_low[:4] in [
+ "krog",
+ "krow",
+ "pcow",
+ "pcog",
+ "pcwg",
+ ]:
+ snu = 1
+ hyst = False
+ if q0_low[-1] == "h":
+ hyst = True
+ q0_low = tokens[0][:-1]
+ if len(q0_low) == 3:
+ what = q0_low[:3]
+ elif q0_low in ["krow", "krog", "pcow", "pcog", "pcwg"]:
+ what = q0_low[:4]
+ elif q0_low[:3] in ["krw", "krg"]:
+ what = q0_low[:3]
+ snu = int(tokens[0][3:])
+ else:
+ what = q0_low[:4]
+ snu = int(tokens[0][4:])
+ if not os.path.isfile(f"{case}.INIT"):
+ plopm_error(
+ f"Missing {cli_error_value(f'{case}.INIT')}, required by "
+ f"{cli_info_value(f'-v {q0_low}')}."
+ )
+ init = OpmFile(f"{case}.INIT")
+ tabdim = init["TABDIMS"]
+ table = np.array(init["TAB"])
+ nswe = tabdim[24]
+ nsnum = tabdim[25]
+ vunit = ""
+ if what == "krg":
+ tunit = "s$_w$ [-]" if use_sw else "s$_g$ [-]"
+ sht = tabdim[23] - 1
+ base = sht + (snu - 1) * nswe
+ time = table[base : base + nswe]
+ time = time[time <= 1.0]
+ count_v = len(time)
+ var = table[
+ sht + nswe * nsnum + (snu - 1) * nswe : sht + nswe * nsnum + snu * nswe
+ ][:count_v]
+ if hyst:
+ base2 = sht + (nsnum // 2 + snu - 1) * nswe
+ timeh = table[base2 : base2 + nswe]
+ timeh = timeh[timeh <= 1.0]
+ count_v = len(timeh)
+ var = np.append(
+ var,
+ np.flip(
+ table[
+ sht
+ + nswe * nsnum
+ + (nsnum // 2 + snu - 1) * nswe : sht
+ + nswe * nsnum
+ + (nsnum // 2 + snu) * nswe
+ ][:count_v]
+ ),
+ )
+ time = np.append(time, np.flip(timeh))
+ if use_sw:
+ time = 1.0 - time
+ elif what == "krow":
+ nswe = tabdim[21]
+ tunit = "s$_w$ [-]"
+ sht = tabdim[26] - 1
+ base = sht + (snu - 1) * nswe
+ time = table[base : base + nswe]
+ time = time[time <= 1.0]
+ count_v = len(time)
+ if tabdim[22] == 2:
+ sht += nswe
+ var = np.flip(
+ table[
+ sht
+ + nswe * nsnum
+ + (snu - 1) * nswe : sht
+ + nswe * nsnum
+ + snu * nswe
+ ][:count_v]
+ )
+ elif what == "krw":
+ nswe = tabdim[21]
+ tunit = "s$_w$ [-]"
+ sht = tabdim[20] - 1
+ base = sht + (snu - 1) * nswe
+ time = table[base : base + nswe]
+ time = time[time <= 1.0]
+ count_v = len(time)
+ var = table[
+ sht + nswe * nsnum + (snu - 1) * nswe : sht + nswe * nsnum + snu * nswe
+ ][:count_v]
+ if hyst:
+ base2 = sht + (nsnum // 2 + snu - 1) * nswe
+ timeh = table[base2 : base2 + nswe]
+ timeh = timeh[timeh <= 1.0]
+ count_v = len(timeh)
+ var = np.append(
+ np.flip(var),
+ table[
+ sht
+ + nswe * nsnum
+ + (nsnum // 2 + snu - 1) * nswe : sht
+ + nswe * nsnum
+ + (nsnum // 2 + snu) * nswe
+ ][:count_v],
+ )
+ time = np.append(np.flip(time), timeh)
+ elif what == "pcow":
+ nswe = tabdim[21]
+ tunit = "s$_w$ [-]"
+ sht = tabdim[20] - 1
+ base = sht + (snu - 1) * nswe
+ time = table[base : base + nswe]
+ time = time[time <= 1.0]
+ count_v = len(time)
+ var = table[
+ sht
+ + 2 * nswe * nsnum
+ + (snu - 1) * nswe : sht
+ + 2 * nswe * nsnum
+ + snu * nswe
+ ][:count_v]
+ else:
+ tunit = "s$_g$ [-]"
+ sht = tabdim[23] - 1
+ base = sht + (snu - 1) * nswe
+ time = table[base : base + nswe]
+ time = time[time <= 1.0]
+ count_v = len(time)
+ var = table[
+ sht
+ + 2 * nswe * nsnum
+ + (snu - 1) * nswe : sht
+ + 2 * nswe * nsnum
+ + snu * nswe
+ ][:count_v]
+ elif values[:6] == "pcfact" or values[:8] == "permfact":
+ cap = 6 if values[:6] == "pcfact" else 8
+ snu = int(tokens[0][cap:]) if not values in ["pcfact", "permfact"] else 1
+ tmp0 = []
+ tmp2 = []
+ found = False
+ vec = tokens[0].upper()[:cap]
+ file_name = _find_keyword(case, vec)
+ count = 0
+ with open(file_name, "r", encoding="utf8") as file:
+ for row in csv.reader(file, delimiter=" "):
+ if len(row) > 0:
+ if row[0] == vec:
+ found = True
+ if count == snu:
+ break
+ if (
+ len(row) > 1
+ and row[0].strip() != "--"
+ and found
+ and count == snu - 1
+ ):
+ tmp0.append(float(row[0]))
+ tmp2.append(float(row[1]))
+ if len(row) > 2 and row[2].strip() == "/":
+ break
+ if (
+ found
+ and row[0] == "/"
+ or len(row) > 2
+ and row[2].strip() == "/"
+ ):
+ count += 1
+ if not tmp2:
+ plopm_error(f"not found {cli_error_value(f'-v {tokens[0]}')}.")
+ var = np.array(tmp2)
+ time = np.array(tmp0)
+ else:
+ summary = OpmSummary(f"{case}.SMSPEC")
+ key = tokens[0].upper()
+ keys = summary.keys()
+ if tokens[0] in cfg.summary_mass:
+ var = summary[key[:-1]]
+ elif key in summary:
+ var = summary[key]
+ else:
+ plopm_error(f"no {cli_error_value(f'-v {tokens[0]}')} found.")
+ if len(tokens) > 1:
+ ops = tokens[1::2]
+ for index, val in enumerate(tokens[2::2]):
+ if val.upper() in keys:
+ other = summary[val.upper()]
+ else:
+ other = float(val)
+ var = _apply_operator(var, other, ops[index])
+ if tunit == "Dates":
+ smsp_dates = 86400 * summary["TIME"]
+ time = np.array(
+ [
+ summary.start_date + datetime.timedelta(seconds=float(sec))
+ for sec in smsp_dates
+ ]
+ )
+ else:
+ time = summary["TIME"] * tskl
+ if tokens[0] in ["fgip", "fgit"]:
+ vunit = " [sm$^3$]"
+ elif tokens[0] in cfg.summary_mass:
+ var *= GAS_DEN_REF
+ vunit = mass_unit(qskl)
+ elif tokens[0] in ["time"]:
+ vunit = " [d]"
+ return time, var * qskl, tunit, vunit
+
+
+
+def _find_keyword(case: str, vec: str) -> str:
+ """Find the deck file containing an OPM keyword.
+
+ Parameters
+ ----------
+ case : str
+ Simulation-case stem.
+ vec : str
+ OPM keyword to locate.
+
+ Returns
+ -------
+ str
+ DATA or included file containing the keyword.
+
+ """
+ include = False
+ path = ""
+ parts = case.split("/")
+ if len(parts) > 1:
+ path = "/".join(parts[:-1]) + "/"
+ case_file = case + ".DATA"
+ includes = []
+ with open(case_file, "r", encoding="utf8") as file:
+ for row in csv.reader(file):
+ if not row:
+ continue
+ val = row[0]
+ if val == vec:
+ return case_file
+ if val == "INCLUDE":
+ include = True
+ continue
+ if include:
+ name = val.split("/")[0].strip(" ")
+ if "'" in name:
+ name = name[1:-1]
+ full = path + name
+ if os.path.isfile(full):
+ includes.append(full)
+ include = False
+ for include_file in includes:
+ with open(include_file, "r", encoding="utf8") as file:
+ for row in csv.reader(file):
+ if not row:
+ continue
+ if row[0] == vec:
+ return include_file
+ files = case_file
+ if len(includes) > 1:
+ if len(includes) == 1:
+ files += f" and {includes[0]}"
+ else:
+ files += ", "
+ files += ", ".join(includes[:-1])
+ files += f" and {includes[-1]}"
+ plopm_error(f"not found keyword {cli_error_value(f'-v {vec}')} " f"in {files}.")
+
+
+def _apply_operator(
+ var: NDArray[np.float64], other: NDArray[np.float64], op: str
+) -> NDArray[np.float64]:
+ """Apply an arithmetic or comparison operator.
+
+ Parameters
+ ----------
+ var, other : np.ndarray
+ Left- and right-hand values.
+ op : str
+ Arithmetic or comparison operator.
+
+ Returns
+ -------
+ np.ndarray
+ Operation result. Failed comparisons are NaN.
+
+ """
+ if op == "+":
+ return var + other
+ if op == "-":
+ return var - other
+ if op == "*":
+ return var * other
+ if op == "/":
+ return var / other
+ mask = ~np.isnan(var)
+ qmask = ~np.isnan(other)
+ mask = mask & qmask
+ if op == "==":
+ var[mask] = np.where(var[mask] == other[mask], 1.0, np.nan)
+ elif op == ">=":
+ var[mask] = np.where(var[mask] >= other[mask], 1.0, np.nan)
+ elif op == "<=":
+ var[mask] = np.where(var[mask] <= other[mask], 1.0, np.nan)
+ elif op == "<":
+ var[mask] = np.where(var[mask] < other[mask], 1.0, np.nan)
+ elif op == ">":
+ var[mask] = np.where(var[mask] > other[mask], 1.0, np.nan)
+ elif op == "!=":
+ var[mask] = np.where(var[mask] != other[mask], 1.0, np.nan)
+ else:
+ plopm_error(f"unknow operation {cli_error_value(f'-v {op}')}.")
+ return var
+
+
+
+[docs]
+def time_unit(times: str) -> tuple[float, str]:
+ """Get the conversion and label for a time unit.
+
+ Parameters
+ ----------
+ times : str
+ Time-unit code or ``"dates"``.
+
+ Returns
+ -------
+ tuple[float, str]
+ Factor converting OPM days and the axis label.
+
+ """
+ if times == "s":
+ return 86400.0, "Time [seconds]"
+ if times == "m":
+ return 1440.0, "Time [minutes]"
+ if times == "h":
+ return 24.0, "Time [hours]"
+ if times == "d":
+ return 1.0, "Time [days]"
+ if times == "w":
+ return 0.14285714285714285, "Time [weeks]"
+ if times == "y":
+ return 0.002737909255898758, "Time [years]"
+ if times == "dates":
+ return 1, "Dates"
+ return 86400.0, "Time [seconds]"
+
+
+
+
+[docs]
+def read_csv_grid(
+ cfg: PlopmConfig, deck: str, n: int
+) -> tuple[NDArray, NDArray, int, int, str, str]:
+ """Read coordinate meshes from a regular CSV grid.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ CSV column and animation configuration.
+ deck : str
+ CSV path without the extension.
+ n : int
+ Map index.
+
+ Returns
+ -------
+ tuple
+ Coordinate meshes, dimensions, and axis names.
+
+ """
+ if cfg.gif:
+ file_name = deck.replace("PLOPM", str(cfg.restart[0]))
+ else:
+ file_name = deck
+ csvv = np.genfromtxt(f"{file_name}.csv", delimiter=",", skip_header=1)
+ col_x = cfg.csv_columns[n][0] - 1
+ col_y = cfg.csv_columns[n][1] - 1
+ x0 = csvv[0, col_x]
+ x1 = csvv[-1, col_x]
+ y0 = csvv[0, col_y]
+ y1 = csvv[-1, col_y]
+ x = x1 + x0
+ y = y1 + y0
+ mx = round(x / (2.0 * x0))
+ my = round(y / (2.0 * y0))
+ xname = "x"
+ yname = "y"
+ xmx = np.linspace(0, x, mx + 1)
+ ymy = np.linspace(0, y, my + 1)
+ return xmx[None, :], ymy[::-1][:, None], mx, my, xname, yname
+
+
+
+def _apply_filter(porvs: NDArray, other: NDArray, op: str, value: float) -> NDArray:
+ """Apply a comparison filter to pore-volume values.
+
+ Parameters
+ ----------
+ porvs : np.ndarray
+ Pore-volume values.
+ other : np.ndarray
+ Values tested by the filter.
+ op : str
+ Comparison operator.
+ value : float
+ Comparison threshold.
+
+ Returns
+ -------
+ np.ndarray
+ Pore volume where the condition is true and zero elsewhere.
+
+ """
+ if op == "==":
+ mask = other == value
+ elif op == ">=":
+ mask = other >= value
+ elif op == "<=":
+ mask = other <= value
+ elif op == "<":
+ mask = other < value
+ elif op == ">":
+ mask = other > value
+ elif op == "!=":
+ mask = other != value
+ else:
+ plopm_error(f"unknow filter operation {cli_error_value(f'-flt {op}')}.")
+ return np.where(mask, porvs, 0)
+
+
+
+[docs]
+def get_unit(name: str) -> str:
+ """Get the display unit for a variable.
+
+ Parameters
+ ----------
+ name : str
+ Variable name.
+
+ Returns
+ -------
+ str
+ Matplotlib-formatted unit label.
+
+ """
+ name_low = name.lower()
+ if name_low in {"disperc", "depth", "dx", "dy", "dz"}:
+ return " [m]"
+ if name_low in {"porv", "fgip", "fgit"}:
+ return r" [sm$^3$]"
+ if name_low in {"permx", "permy", "permz"}:
+ return " [mD]"
+ if name_low in {"tranx", "trany", "tranz"}:
+ return " [cP rm$^3$/ (day bar)]"
+ if name_low in {"pressure", "rpr", "fpr", "fprr", "wbhp"}:
+ return " [bar]"
+ return " [-]"
+
+
+
+
+[docs]
+def read_quantity(
+ deck: str,
+ data: SimData,
+ name: str,
+ step: int,
+ scale: float,
+ mass: list[str],
+ mass_all: list[str],
+ caprock: list[str],
+ stress: float,
+ filters: str,
+ isgif: bool,
+ vmin: str,
+ vmax: str,
+ cvs: list,
+) -> tuple[str, NDArray]:
+ """Read and transform one spatial or VTK quantity.
+
+ Parameters
+ ----------
+ deck : str
+ Simulation-case stem or CSV path.
+ data : SimData
+ Loaded simulation data.
+ name : str
+ Variable name or expression.
+ step : int
+ Restart report step.
+ scale : float
+ Scale factor applied to derived values.
+ mass, mass_all, caprock : list[str]
+ Supported derived-variable groups.
+ stress : float
+ Stress coefficient for caprock quantities.
+ filters : str
+ Property-filter expression.
+ isgif : bool
+ Whether the CSV path contains a restart placeholder.
+ vmin, vmax : str
+ Optional value thresholds.
+ cvs : list
+ CSV input and column settings.
+
+ Returns
+ -------
+ tuple[str, np.ndarray]
+ Unit label and quantity values.
+
+ """
+ names = name.split(" ")
+ unit = get_unit(name)
+ name0_low = names[0]
+ name0 = name0_low.upper()
+ if cvs[0]:
+ if isgif:
+ file_name = deck.replace("PLOPM", str(step))
+ else:
+ file_name = deck
+ csvv = np.genfromtxt(f"{file_name}.csv", delimiter=",", skip_header=1)
+ col = cvs[2] - 1
+ values = csvv[:, col]
+ else:
+ if data.init.count(name0):
+ values = np.array(data.init[name0], dtype=float)
+ if name0_low == "porv":
+ values = data.active_pv
+ elif name0_low in ["wells", "faults", "grid"]:
+ values = np.zeros_like(data.init["SATNUM"])
+ elif name0_low in ["index_i", "index_j", "index_k"]:
+ values = np.array(
+ _grid_indices(name0_low, data.nx, data.ny, data.nz), dtype=float
+ )
+ values = values[data.porv > 0]
+ elif data.unrst.count(name0, step):
+ values = data.unrst[name0, step]
+ if data.unrst.count("RPORV", step):
+ if filters:
+ porv0 = np.array(data.init["PORV"])
+ mask = porv0 > 0
+ base_rporv = np.array(data.unrst["RPORV", step])
+ for value in filters.split("&"):
+ filte = value.strip().split(" ")
+ key = filte[0].upper()
+ if data.init.count(key):
+ q1 = np.array(data.init[key])
+ elif data.unrst.count(key, step):
+ q1 = np.array(data.unrst[key, step])
+ else:
+ plopm_error(
+ f"unknow filter quantity {cli_error_value(f'-flt {key}')}."
+ )
+ base_rporv = _apply_filter(
+ base_rporv, q1, filte[1], float(filte[2])
+ )
+ data.porv[mask] = base_rporv
+ else:
+ data.porv[data.porv > 0] = np.array(data.unrst["RPORV", step])
+ elif name0_low in mass_all:
+ values = _get_mass(data, name0_low, step) * scale
+ if name0_low in mass:
+ unit = mass_unit(scale)
+ elif name0_low in caprock:
+ values, unit = _get_caprock(data, name0_low, step, stress)
+ elif name0_low in ["swat", "soil", "sgas"]:
+ values = _get_saturation(data.unrst, name0_low, step) * scale
+ else:
+ plopm_error(f"not found {cli_error_value(f'-v {name0}')}.")
+ if len(names) > 1:
+ ops = names[1::2]
+ for j, val in enumerate(names[2::2]):
+ if val[0].isdigit() and not val[-1].isdigit():
+ q1 = data.unrst[val[1:].upper(), int(val[0])]
+ elif val[0].isdigit() and val[-1].isdigit():
+ q1 = np.full_like(values, float(val))
+ elif data.init.count(val.upper()):
+ q1 = np.array(data.init[val.upper()])
+ if val.upper() == "PORV":
+ q1 = q1[data.porv > 0]
+ elif val in ["index_i", "index_j", "index_k"]:
+ q1 = np.array(
+ _grid_indices(val, data.nx, data.ny, data.nz),
+ dtype=float,
+ )
+ q1 = q1[data.porv > 0]
+ elif data.unrst.count(val.upper(), step):
+ q1 = data.unrst[val.upper(), step]
+ elif val in mass_all:
+ q1 = _get_mass(data, val, step) * scale
+ elif val in caprock:
+ q1, unit = _get_caprock(data, val, step, stress)
+ else:
+ plopm_error(f"not found {cli_error_value(f'-v {val}')}.")
+ values = _apply_operator(values, q1, ops[j])
+ if vmin:
+ values = np.asarray(values)
+ values[values < float(vmin)] = np.nan
+ if vmax:
+ values = np.asarray(values)
+ values[values > float(vmax)] = np.nan
+ return unit, values
+
+
+
+def _get_saturation(unrst: OpmRestart, name: str, step: int) -> NDArray:
+ """Derive a missing phase saturation.
+
+ Parameters
+ ----------
+ unrst : OpmRestart
+ UNRST reader.
+ name : {"soil", "swat", "sgas"}
+ Saturation to derive.
+ step : int
+ Restart report step.
+
+ Returns
+ -------
+ np.ndarray
+ Requested phase saturation.
+
+ """
+ if unrst.count("SOIL", step):
+ soil = np.array(unrst["SOIL", step])
+ else:
+ soil = np.array(0)
+ if unrst.count("SGAS", step):
+ sgas = np.array(unrst["SGAS", step])
+ else:
+ sgas = np.array(0)
+ if unrst.count("SWAT", step):
+ swat = np.array(unrst["SWAT", step])
+ else:
+ swat = np.array(0)
+ if name == "soil":
+ return 1 - sgas - swat
+ if name == "swat":
+ return 1 - sgas - soil
+ return 1 - soil - swat
+
+
+def _get_mass(data: SimData, name: str, step: int) -> NDArray:
+ """Compute component masses and mass fractions.
+
+ Parameters
+ ----------
+ data : SimData
+ Loaded restart data and pore volume.
+ name : str
+ Requested derived variable.
+ step : int
+ Restart report step.
+
+ Returns
+ -------
+ np.ndarray
+ Requested component quantity.
+
+ """
+ sgas = np.array(data.unrst["SGAS", step])
+ rhog = np.array(data.unrst["GAS_DEN", step])
+ rhow = np.array(data.unrst["WAT_DEN", step])
+ if data.unrst.count("RSW", step):
+ rsw = np.array(data.unrst["RSW", step])
+ else:
+ rsw = np.zeros_like(sgas)
+ if data.unrst.count("RVW", step):
+ rvw = np.array(data.unrst["RVW", step])
+ else:
+ rvw = np.zeros_like(sgas)
+ if data.unrst.count("RPORV", step):
+ rpv = np.array(data.unrst["RPORV", step])
+ else:
+ rpv = data.active_pv
+ denom_l = rsw + WAT_DEN_REF / GAS_DEN_REF
+ denom_g = rvw + GAS_DEN_REF / WAT_DEN_REF
+ x_l_co2 = np.zeros_like(rsw)
+ x_g_h2o = np.zeros_like(rvw)
+ mask_l = denom_l != 0
+ mask_g = denom_g != 0
+ x_l_co2[mask_l] = rsw[mask_l] / denom_l[mask_l]
+ x_g_h2o[mask_g] = rvw[mask_g] / denom_g[mask_g]
+ inv_sgas = 1.0 - sgas
+ inv_xg = 1.0 - x_g_h2o
+ inv_xl = 1.0 - x_l_co2
+ co2_g = inv_xg * sgas * rhog * rpv
+ co2_d = x_l_co2 * inv_sgas * rhow * rpv
+ h2o_l = inv_xl * inv_sgas * rhow * rpv
+ h2o_v = x_g_h2o * sgas * rhog * rpv
+ return _select_mass(name, co2_g, co2_d, h2o_l, h2o_v, x_l_co2, x_g_h2o)
+
+
+def _select_mass(
+ name: str,
+ co2_g: NDArray,
+ co2_d: NDArray,
+ h2o_l: NDArray,
+ h2o_v: NDArray,
+ x_l_co2: NDArray,
+ x_g_h2o: NDArray,
+) -> NDArray:
+ """Select a mass or mass-fraction result by name.
+
+ Parameters
+ ----------
+ name : str
+ Requested derived variable.
+ co2_g, co2_d : np.ndarray
+ Free and dissolved CO2 masses.
+ h2o_l, h2o_v : np.ndarray
+ Liquid and vapor water masses.
+ x_l_co2, x_g_h2o : np.ndarray
+ CO2-in-liquid and water-in-gas mass fractions.
+
+ Returns
+ -------
+ np.ndarray
+ Selected mass or mass fraction.
+
+ """
+ if name == "gasm":
+ return co2_g
+ if name == "dism":
+ return co2_d
+ if name == "liqm":
+ return h2o_l
+ if name == "vapm":
+ return h2o_v
+ if name == "h2om":
+ return h2o_v + h2o_l
+ if name == "xco2l":
+ return x_l_co2
+ if name == "xh2ov":
+ return x_g_h2o
+ if name == "xco2v":
+ return 1 - x_g_h2o
+ if name == "xh2ol":
+ return 1 - x_l_co2
+ return co2_g + co2_d
+
+
+def _get_caprock(
+ data: SimData, name: str, step: int, stress: float
+) -> tuple[NDArray, str]:
+ """Compute a caprock-integrity quantity.
+
+ Parameters
+ ----------
+ data : SimData
+ Loaded static and restart properties.
+ name : str
+ Requested caprock variable.
+ step : int
+ Restart report step.
+ stress : float
+ Vertical stress coefficient.
+
+ Returns
+ -------
+ tuple[np.ndarray, str]
+ Computed values and unit label.
+
+ """
+ init = data.init
+ unrst = data.unrst
+ dz = np.array(init["DZ", 0])
+ depth = np.array(init["DEPTH", 0])
+ dz_half = 0.5 * dz
+ dz_corr = 0.5 * dz
+ if unrst.count("WAT_DEN", 0) and unrst.count("WAT_DEN", step):
+ den0 = np.array(unrst["WAT_DEN", 0])
+ den1 = np.array(unrst["WAT_DEN", step])
+ else:
+ den0 = np.array(1000.0)
+ den1 = np.array(1000.0)
+ fac = 9.81 / 1e5
+ pz_c0 = fac * dz_corr * den0
+ pz_c1 = fac * dz_corr * den1
+ pressure0 = np.array(unrst["PRESSURE", 0])
+ pressure1 = np.array(unrst["PRESSURE", step])
+ limipres = stress * (depth - dz_half)
+ overpres = limipres - (pressure1 - pz_c1)
+ limipres -= pressure0 - pz_c0
+ objepres = np.zeros_like(overpres)
+ mask = limipres != 0
+ objepres[mask] = overpres[mask] / limipres[mask]
+ if name == "limipres":
+ return limipres, " [bar]"
+ if name == "overpres":
+ return -overpres, " [bar]"
+ return objepres, " [-]"
+
+
+
+[docs]
+def get_wells(cfg: PlopmConfig, n: int) -> tuple[list, list]:
+ """Read wells intersecting the selected slice.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Case and slice configuration.
+ n : int
+ Case index.
+
+ Returns
+ -------
+ tuple[list, list[str]]
+ Completion intervals grouped by well and the well names.
+
+ """
+ wells: list[list[list[int]]] = []
+ lwells: list[str] = []
+ well_map = {}
+ haswells = False
+ sources = False
+ with open(f"{cfg.cases[0][n]}.DATA", "r", encoding="utf8") as file:
+ for row in csv.reader(file):
+ if not row:
+ continue
+ tokens = row[0].split()
+ if not tokens:
+ continue
+ key = tokens[0]
+ if key == "COMPDAT":
+ haswells = True
+ continue
+ if key == "SOURCE":
+ sources = True
+ continue
+ if key == "/":
+ haswells = False
+ sources = False
+ continue
+ if key.startswith("--"):
+ continue
+ if haswells:
+ if len(tokens) < 5:
+ continue
+ wname = tokens[0]
+ if wname not in well_map:
+ well_map[wname] = len(lwells)
+ lwells.append(wname)
+ wells.append([])
+ idx = well_map[wname]
+ wells[idx].append(
+ [
+ int(tokens[1]) - 1,
+ int(tokens[2]) - 1,
+ int(tokens[3]) - 1,
+ int(tokens[4]) - 1,
+ ]
+ )
+ elif sources:
+ if len(tokens) < 3:
+ continue
+ wname = tokens[0]
+ if wname not in well_map:
+ well_map[wname] = len(lwells)
+ lwells.append(wname)
+ wells.append([])
+ idx = well_map[wname]
+ wells[idx].append(
+ [
+ int(tokens[0]) - 1,
+ int(tokens[1]) - 1,
+ int(tokens[2]) - 1,
+ int(tokens[2]) - 1,
+ ]
+ )
+ if not cfg.global_range:
+ sld_x = cfg.slice[n][0]
+ sld_y = cfg.slice[n][1]
+ sld_z = cfg.slice[n][2]
+ whow = cfg.slice_mode
+ for i, wells_list in enumerate(wells):
+ for j, well in enumerate(wells_list):
+ if not well:
+ continue
+ keep = True
+ if sld_x[0] > -1:
+ val = well[0]
+ if whow == "min":
+ keep = sld_x[0] <= val < sld_x[1]
+ else:
+ keep = val == sld_x[0]
+ elif sld_y[0] > -1:
+ val = well[1]
+ if whow == "min":
+ keep = sld_y[0] <= val < sld_y[1]
+ else:
+ keep = val == sld_y[0]
+ else:
+ z0, z1 = well[2], well[3]
+ if whow == "min":
+ keep = not (sld_z[1] < z0 or sld_z[0] > z1)
+ else:
+ keep = sld_z[0] >= z0 and sld_z[0] <= z1
+ if not keep:
+ wells[i][j] = []
+ return wells, lwells
+
+
+
+
+[docs]
+def get_faults(cfg: PlopmConfig, n: int) -> tuple[list, list]:
+ """Read faults intersecting the selected slice.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Case and slice configuration.
+ n : int
+ Case index.
+
+ Returns
+ -------
+ tuple[list, list[str]]
+ Grid segments grouped by fault and the fault names.
+
+ """
+ faults: list[list[list[int]]] = []
+ lfaults: list[str] = []
+ fault_map = {}
+ hasfaults = False
+ with open(f"{cfg.cases[0][n]}.DATA", "r", encoding="utf8") as file:
+ for row in csv.reader(file):
+ if not row:
+ continue
+ tokens = row[0].split()
+ if not tokens:
+ continue
+ key = tokens[0]
+ if key == "FAULTS":
+ hasfaults = True
+ continue
+ if hasfaults:
+ if key.startswith("--"):
+ continue
+ if "/" in key:
+ break
+ if len(tokens) < 7:
+ continue
+ fname = key
+ if fname not in fault_map:
+ fault_map[fname] = len(lfaults)
+ lfaults.append(fname)
+ faults.append([])
+ idx = fault_map[fname]
+ faults[idx].append(
+ [
+ int(tokens[1]) - 1,
+ int(tokens[3]) - 1,
+ int(tokens[5]) - 1,
+ int(tokens[6]) - 1,
+ ]
+ )
+ if not cfg.global_range:
+ sld_x = cfg.slice[n][0]
+ sld_y = cfg.slice[n][1]
+ sld_z = cfg.slice[n][2]
+ whow = cfg.slice_mode
+ for i, flist in enumerate(faults):
+ for j, fault in enumerate(flist):
+ if not fault:
+ continue
+ keep = True
+ if sld_x[0] > -1:
+ val = fault[0]
+ if whow == "min":
+ keep = sld_x[0] <= val < sld_x[1]
+ else:
+ keep = val == sld_x[0]
+ elif sld_y[0] > -1:
+ val = fault[1]
+ if whow == "min":
+ keep = sld_y[0] <= val < sld_y[1]
+ else:
+ keep = val == sld_y[0]
+ else:
+ z0, z1 = fault[2], fault[3]
+ if whow == "min":
+ keep = not (sld_z[1] < z0 or sld_z[0] > z1)
+ else:
+ keep = sld_z[0] >= z0 and sld_z[0] <= z1
+ if not keep:
+ faults[i][j] = []
+ return faults, lfaults
+
+
+
+
+
+
+# SPDX-FileCopyrightText: 2026 NORCE Research AS
+# SPDX-License-Identifier: GPL-3.0
+
+"""Format plopm help text and command-line messages.
+
+The module hides deprecated aliases from ``--help``, reports their replacements,
+and applies ANSI colors only when supported by the selected output stream.
+"""
+
+import argparse
+import os
+import sys
+from collections.abc import Sequence
+from typing import NoReturn
+
+DEPRECATED_OPTION_ALIASES = {
+ # Input and output
+ "-csv": "-cc",
+ "--csv": "--csv-columns",
+ "-p": "-fp",
+ "--path": "--flow-path",
+ "--mode": "--format",
+ "--output": "--output-dir",
+ "-save": "-fn",
+ "--save": "--filename",
+ # Spatial and temporal selection
+ "--slide": "--slice",
+ "-tunits": "-tu",
+ "--tunits": "--time-units",
+ "-distance": "-dist",
+ # Filtering, masking, and thresholds
+ "-filter": "-flt",
+ "--vmin": "--min-threshold",
+ "--vmax": "--max-threshold",
+ "-mask": "-mv",
+ "--mask": "--mask-variable",
+ "-maskthr": "-mt",
+ "--maskthr": "--mask-threshold",
+ # Computation and data transformation
+ "-how": "-agg",
+ "--how": "--aggregation",
+ "-a": "-sf",
+ "--adjust": "--scale-factor",
+ "-diff": "-di",
+ "--diff": "--difference-input",
+ "-stress": "-sc",
+ "--stress": "--stress-coefficient",
+ "-dual": "-dg",
+ "--dual": "--dual-grid",
+ # Plot types and statistical representation
+ "-histogram": "-hist",
+ "-ensemble": "-ens",
+ "-bandprop": "-fb",
+ "--bandprop": "--fill-between-style",
+ "-step": "-sp",
+ "--step": "--step-plot",
+ # Figure and subplot layout
+ "-d": "-fs",
+ "--dimensions": "--figsize",
+ "-subfigs": "-sg",
+ "--subfigs": "--subplot-grid",
+ "-cbsfax": "-cbp",
+ "--cbsfax": "--colorbar-position",
+ "-delax": "-rdl",
+ "--delax": "--remove-duplicate-labels",
+ # Titles, labels, and legends
+ "-suptitle": "-st",
+ "-xlabel": "-xl",
+ "-ylabel": "-yl",
+ "-clabel": "-cbl",
+ "--clabel": "--colorbar-label",
+ "-labels": "-llb",
+ "--labels": "--legend-labels",
+ "-loc": "-ll",
+ "--loc": "--legend-location",
+ "-remove": "-hide",
+ "--remove": "--hide-map-elements",
+ # Axes, coordinates, and formatting
+ "-xunits": "-xu",
+ "-yunits": "-yu",
+ "-z": "-asp",
+ "--scale": "--equal-aspect",
+ "-rotate": "-rot",
+ "--rotate": "--rotation",
+ "-translate": "-tr",
+ "--translate": "--translation",
+ "-xformat": "-xf",
+ "-yformat": "-yf",
+ "-xlnum": "-xnt",
+ "--xlnum": "--xtick-count",
+ "-ylnum": "-ynt",
+ "--ylnum": "--ytick-count",
+ # Color scales and styling
+ "-b": "-cl",
+ "--bounds": "--clim",
+ "-log": "-clog",
+ "--log": "--color-log",
+ "-clogthks": "-clt",
+ "--clogthks": "--color-log-ticks",
+ "-global": "-gr",
+ "--global": "--global-range",
+ "-cformat": "-cbf",
+ "--cformat": "--colorbar-format",
+ "-cnum": "-cbn",
+ "--cnum": "--colorbar-tick-count",
+ "-cticks": "-cbt",
+ "--cticks": "--colorbar-ticks",
+ "--lw": "--linewidth",
+ "-e": "-ls",
+ "-axgrid": "-ag",
+ "--axgrid": "--axis-grid",
+ "-facecolor": "-fc",
+ "-ncolor": "-ic",
+ "--ncolor": "--inactive-color",
+ "-grid": "-ge",
+ "--grid": "--grid-edges",
+ "-f": "-fz",
+ "--size": "--fontsize",
+ # VTK output
+ "-vtkformat": "-vf",
+ "--vtkformat": "--vtk-format",
+ "-vtknames": "-vn",
+ "--vtknames": "--vtk-names",
+ # GIF output
+ "-interval": "-gi",
+ "--interval": "--gif-interval",
+ "-loop": "-gl",
+ "--loop": "--gif-loop",
+ # Information and diagnostics
+ "-printv": "-lv",
+ "--printv": "--list-variables",
+}
+ANSI_BOLD_RED = "1;31"
+ANSI_BOLD_YELLOW = "1;33"
+ANSI_BOLD_GREEN = "1;32"
+ANSI_BOLD_BLUE = "1;34"
+ANSI_BOLD_MAGENTA = "1;35"
+ANSI_YELLOW = "1;33"
+ANSI_GREEN = "1;32"
+ANSI_CYAN = "36"
+ANSI_RED = "31"
+ANSI_BLUE = "1;34"
+
+
+
+[docs]
+class PlopmHelpFormatter(argparse.ArgumentDefaultsHelpFormatter):
+ """Argparse formatter that hides deprecated option aliases.
+
+ Current options retain the standard
+ :class:`argparse.ArgumentDefaultsHelpFormatter` layout and default values.
+
+ """
+
+ def _format_action_invocation(
+ self,
+ action: argparse.Action,
+ ) -> str:
+ """Format one argparse action using current option names.
+
+ Parameters
+ ----------
+ action : argparse.Action
+ Action whose option invocation is displayed in CLI help.
+
+ Returns
+ -------
+ str
+ Formatted invocation with deprecated aliases omitted.
+
+ """
+ if not action.option_strings:
+ return super()._format_action_invocation(action)
+
+ original_options = action.option_strings
+ visible_options = [
+ option
+ for option in original_options
+ if option not in DEPRECATED_OPTION_ALIASES
+ ]
+
+ if not visible_options:
+ return super()._format_action_invocation(action)
+
+ try:
+ action.option_strings = visible_options
+ return super()._format_action_invocation(action)
+ finally:
+ action.option_strings = original_options
+
+
+
+
+[docs]
+def warn_deprecated_options(argv: Sequence[str]) -> None:
+ """Warn once for each deprecated option in an argument list.
+
+ Parameters
+ ----------
+ argv : Sequence[str]
+ Command-line arguments, excluding or including the executable name.
+
+ """
+ reported: set[str] = set()
+
+ for argument in argv:
+ option = argument.partition("=")[0]
+
+ if option not in DEPRECATED_OPTION_ALIASES or option in reported:
+ continue
+
+ replacement = DEPRECATED_OPTION_ALIASES[option]
+
+ plopm_warning(
+ f"option {cli_deprecated_value(option)} is deprecated and will be "
+ f"removed in the next release; use {cli_current_value(replacement)} "
+ "instead"
+ )
+ reported.add(option)
+
+
+
+def _supports_color(stream: object = sys.stderr) -> bool:
+ """Check whether an output stream supports ANSI colors.
+
+ Parameters
+ ----------
+ stream : object, default: sys.stderr
+ Output stream to inspect.
+
+ Returns
+ -------
+ bool
+ ``True`` for an interactive stream unless colors are disabled by
+ ``NO_COLOR`` or ``TERM=dumb``.
+
+ """
+ return (
+ hasattr(stream, "isatty")
+ and stream.isatty()
+ and os.environ.get("NO_COLOR") is None
+ and os.environ.get("TERM") != "dumb"
+ )
+
+
+def _colorize(
+ text: str,
+ code: str,
+ stream: object = sys.stderr,
+) -> str:
+ """Wrap text in an ANSI color sequence when supported.
+
+ Parameters
+ ----------
+ text : str
+ Text to format.
+ code : str
+ ANSI Select Graphic Rendition code.
+ stream : object, default: sys.stderr
+ Output stream used to determine color support.
+
+ Returns
+ -------
+ str
+ Colored text, or the original text when colors are unavailable.
+
+ """
+ if not _supports_color(stream):
+ return text
+ return f"\033[{code}m{text}\033[0m"
+
+
+
+[docs]
+def cli_deprecated_value(value: str) -> str:
+ """Format a deprecated CLI option or value.
+
+ Parameters
+ ----------
+ value : str
+ Option or value to display.
+
+ Returns
+ -------
+ str
+ Quoted value with deprecated-option styling when supported.
+
+ """
+ return _colorize(repr(value), ANSI_YELLOW)
+
+
+
+
+[docs]
+def cli_current_value(value: str) -> str:
+ """Format a current CLI option or value.
+
+ Parameters
+ ----------
+ value : str
+ Option or value to display.
+
+ Returns
+ -------
+ str
+ Quoted value with current-option styling when supported.
+
+ """
+ return _colorize(repr(value), ANSI_GREEN)
+
+
+
+
+[docs]
+def cli_error_value(value: str) -> str:
+ """Format an invalid CLI option or value.
+
+ Parameters
+ ----------
+ value : str
+ Option or value to display.
+
+ Returns
+ -------
+ str
+ Quoted value with error styling when supported.
+
+ """
+ return _colorize(repr(value), ANSI_RED)
+
+
+
+
+[docs]
+def cli_info_value(value: str) -> str:
+ """Format an informational CLI option or value.
+
+ Parameters
+ ----------
+ value : str
+ Option or value to display.
+
+ Returns
+ -------
+ str
+ Quoted value with informational styling when supported.
+
+ """
+ return _colorize(repr(value), ANSI_BLUE)
+
+
+
+
+[docs]
+def plopm_error(message: str) -> NoReturn:
+ """Raise a fatal command-line error.
+
+ Parameters
+ ----------
+ message : str
+ Error message displayed after the plopm label.
+
+ Raises
+ ------
+ SystemExit
+ Always raised with the formatted error message.
+
+ """
+ label = _colorize("error", ANSI_BOLD_RED)
+ raise SystemExit(f"{plopm_name()}: {label}: {message}")
+
+
+
+
+[docs]
+def plopm_warning(message: str) -> None:
+ """Display a non-fatal command-line warning.
+
+ Parameters
+ ----------
+ message : str
+ Warning message displayed on standard error.
+
+ """
+ label = _colorize("warning", ANSI_BOLD_YELLOW)
+ print(f"{plopm_name()}: {label}: {message}", file=sys.stderr)
+
+
+
+
+[docs]
+def plopm_info(message: str) -> None:
+ """Display an informational command-line message.
+
+ Parameters
+ ----------
+ message : str
+ Message displayed on standard output.
+
+ """
+ label = _colorize("info", ANSI_BOLD_BLUE, sys.stdout)
+ print(f"{plopm_name()}: {label}: {message}")
+
+
+
+
+[docs]
+def plopm_tip(message: str) -> None:
+ """Display a command-line suggestion.
+
+ Parameters
+ ----------
+ message : str
+ Suggestion displayed on standard output.
+
+ """
+ label = _colorize("tip", ANSI_BOLD_MAGENTA, sys.stdout)
+ print(f"{plopm_name(sys.stdout)}: {label}: {message}")
+
+
+
+
+[docs]
+def plopm_success(output_dir: str, filenames: list[str]) -> None:
+ """Display the generated output location and filenames.
+
+ Parameters
+ ----------
+ output_dir : str
+ Directory containing the generated files.
+ filenames : list[str]
+ Generated filenames.
+
+ """
+ label = _colorize("success", ANSI_BOLD_GREEN, sys.stdout)
+ if not filenames:
+ plopm_error("Unreachable code executed")
+ elif len(filenames) == 1:
+ print(f"{plopm_name()}: {label}: {output_dir}/{filenames[0]}")
+ elif len(filenames) <= 5:
+ print(f"{plopm_name()}: {label}")
+ print(f" Output directory: {output_dir}")
+ print(f" Files: {', '.join(filenames)}")
+ else:
+ print(f"{plopm_name()}: {label}")
+ print(f" Output directory: {output_dir}")
+ print(f" Files ({len(filenames)}):")
+ for filename in filenames:
+ print(f" - {filename}")
+
+
+
+
+[docs]
+def plopm_name(stream: object = sys.stderr) -> str:
+ """Format the plopm program name.
+
+ Parameters
+ ----------
+ stream : object, default: sys.stderr
+ Output stream used to determine color support.
+
+ Returns
+ -------
+ str
+ Program name with gradient colors when supported.
+
+ """
+ characters = [
+ ("p", "36"),
+ ("l", "36"),
+ ("o", "35"),
+ ("p", "36"),
+ ("m", "36"),
+ ]
+ return "".join(
+ _colorize(character, color, stream) for character, color in characters
+ )
+
+
+
+
+
+
+# SPDX-FileCopyrightText: 2024-2026 NORCE Research AS
+# SPDX-License-Identifier: GPL-3.0
+# pylint: disable=W3301,W0123,R0912,R0915,R0914,R1702,W0611,R0913,R0917,C0302,C0115,R0916,E1102
+
+"""Create one-dimensional plots and tabular output from OPM results.
+
+The module reads summary vectors, grid-derived series, and optional CSV data.
+It also supports ensemble statistics, subplot layouts, and PNG or CSV output.
+"""
+
+import os
+import warnings
+
+import matplotlib.pyplot as plt
+import numpy as np
+from matplotlib.axes import Axes
+from matplotlib.figure import Figure
+from numpy.typing import NDArray
+from scipy.interpolate import interp1d
+from scipy.stats import lognorm, norm
+
+from plopm.config.config import PlopmConfig
+from plopm.utils.readers import read_series
+from plopm.utils.terminal import cli_info_value, plopm_info
+
+
+
+[docs]
+def make_plots(cfg: PlopmConfig) -> list[str]:
+ """Create the requested one-dimensional plots and CSV files.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized plotting configuration.
+
+ Returns
+ -------
+ list[str]
+ Names of the generated files.
+
+ """
+ generated_files: list[str] = []
+
+ deckn = _get_deck_name(cfg.cases[0][0])
+ fig, _ = plt.subplots(1, 1)
+ if (
+ cfg.ensemble == 0
+ and not cfg.subplot_grid[0]
+ and len(cfg.cases[0]) < len(cfg.variables)
+ ):
+ cfg.cases[0] = [cfg.cases[0][0]] * len(cfg.variables)
+ if len(cfg.linewidth[0]) < len(cfg.variables):
+ cfg.linewidth[0] = [cfg.linewidth[0][0]] * len(cfg.variables)
+ cfg.linewidth = [cfg.linewidth[0]] * len(cfg.variables)
+ if len(cfg.colors[0]) < len(cfg.variables):
+ cfg.colors[0] = [cfg.colors[0][0]] * len(cfg.variables)
+ cfg.colors = [cfg.colors[0]] * len(cfg.variables)
+ if len(cfg.linestyle[0]) < len(cfg.variables):
+ cfg.linestyle[0] = [cfg.linestyle[0][0]] * len(cfg.variables)
+ cfg.linestyle = [cfg.linestyle[0]] * len(cfg.variables)
+ if cfg.subplot_grid[0]:
+ plt.close()
+ fig, axes = plt.subplots(
+ int(cfg.subplot_grid[0]), int(cfg.subplot_grid[1]), layout="compressed"
+ )
+ for j, quan in enumerate(cfg.variables):
+ k = j
+ if not cfg.subplot_grid[0]:
+ plt.close()
+ fig, axes = plt.subplots(1, 1, layout="compressed")
+ axes = np.array([axes])
+ k = 0
+ axis = axes.flat[k]
+ axis.grid(int(cfg.axis_grid[j]))
+ if cfg.ensemble > 0:
+ tunit, vunit, min_t, max_t, min_v, max_v = _plot_ensemble(cfg, axes)
+ else:
+ ylow = 0 if cfg.ylog[j] == "1" else -np.inf
+ xlow = 0 if cfg.xlog[j] == "1" else -np.inf
+ min_t, max_t, min_v, max_v = 0, 0, 0, 0
+ for i, name in enumerate(cfg.cases[j]):
+ jj = j
+ if len(cfg.variables) == len(cfg.cases[0]) and not cfg.subplot_grid[0]:
+ jj = i
+ quan = cfg.variables[i]
+ time, var, tunit, vunit = read_series(
+ cfg, name, quan, cfg.time_units[jj], float(cfg.scale_factor[jj]), i
+ )
+ label = _get_label(cfg, name, jj, i)
+ if cfg.step_plot:
+ axis.step(
+ time,
+ var,
+ color=cfg.colors[jj][i % len(cfg.colors[jj])],
+ ls=cfg.linestyle[jj][i % len(cfg.linestyle[jj])],
+ label=label,
+ lw=float(cfg.linewidth[jj][i]),
+ )
+ elif cfg.histogram[0]:
+ ij = i + j * len(cfg.cases[j])
+ if (
+ len(cfg.variables) == len(cfg.cases[0])
+ and not cfg.subplot_grid[0]
+ ):
+ ij = i
+ hist = cfg.histogram[ij].split(",")
+ mean = np.nanmean(var)
+ std = np.nanstd(var)
+ plopm_info(
+ f"histogram: {cli_info_value(f'mean={mean:.6E}')}, "
+ f"{cli_info_value(f'std={std:.6E}')}"
+ )
+ if not cfg.legend_labels[0][0]:
+ label += f" (mean={mean:.3E}, std={std:.3E})"
+ counts, bins, _ = axis.hist(
+ var,
+ int(hist[0]),
+ color=cfg.colors[jj][(i + k) % len(cfg.colors[jj])],
+ label=label,
+ )
+ if len(hist) > 1:
+ xnorm = np.linspace(bins[0], bins[-1], 1000)
+ if hist[1] == "norm":
+ norm_pdf = norm.pdf(xnorm, mean, std)
+ norm_max = np.max(norm_pdf)
+ if norm_max > 0:
+ axis.plot(
+ xnorm,
+ np.max(counts) * norm_pdf / norm_max,
+ color=cfg.colors[jj][(i + k) % len(cfg.colors[jj])],
+ )
+ elif hist[1] == "lognorm":
+ if mean > 0:
+ a = 1 + (std / mean) ** 2
+ s = np.sqrt(np.log(a))
+ scale = mean / np.sqrt(a)
+ dist = lognorm(s, 0, scale)
+ dist_pdf = dist.pdf(xnorm)
+ dist_max = np.max(dist_pdf)
+ plopm_info(
+ f"distribution: "
+ f"{cli_info_value(f'lognorm({s:.6E}, 0, {scale:.6E})')}"
+ )
+ if dist_max > 0:
+ axis.plot(
+ xnorm,
+ np.max(counts) * dist_pdf / dist_max,
+ color=cfg.colors[jj][
+ (i + k) % len(cfg.colors[jj])
+ ],
+ )
+ else:
+ axis.plot(
+ time,
+ var,
+ color=cfg.colors[jj][i % len(cfg.colors[jj])],
+ ls=cfg.linestyle[jj][i % len(cfg.linestyle[jj])],
+ label=label,
+ lw=float(cfg.linewidth[jj][i]),
+ )
+ min_t, max_t, min_v, max_v = _update_limits(
+ time, var, tunit, min_t, max_t, min_v, max_v, xlow, ylow, i == 0
+ )
+ axis.set_ylabel(quan + vunit)
+ if not cfg.histogram[0]:
+ if min_v != max_v:
+ axis.set_ylim([min_v, max_v])
+ else:
+ axis.set_ylabel("Histogram of " + quan + vunit)
+ if not cfg.remove_duplicate_labels or k + int(cfg.subplot_grid[1]) >= len(
+ cfg.variables
+ ):
+ axis.set_xlabel(tunit)
+ if cfg.xlabel[0]:
+ axis.set_xlabel(cfg.xlabel[j])
+ if cfg.ylabel[0]:
+ axis.set_ylabel(cfg.ylabel[j])
+ xlabels = np.empty(0)
+ ylabels = np.empty(0)
+ if len(cfg.xlim[0]) > 1:
+ axis.set_xlim([float(cfg.xlim[j][0][1:]), float(cfg.xlim[j][1][:-1])])
+ xlabels = np.linspace(
+ float(cfg.xlim[j][0][1:]),
+ float(cfg.xlim[j][1][:-1]),
+ int(cfg.xtick_count[j]),
+ )
+ elif tunit != "Dates" and not cfg.histogram[0]:
+ if min_v != max_v:
+ axis.set_xlim([min_t, max_t])
+ xlabels = np.linspace(min_t, max_t, int(cfg.xtick_count[j]))
+ if len(cfg.ylim[0]) > 1:
+ axis.set_ylim([float(cfg.ylim[j][0][1:]), float(cfg.ylim[j][1][:-1])])
+ ylabels = np.linspace(
+ float(cfg.ylim[j][0][1:]),
+ float(cfg.ylim[j][1][:-1]),
+ int(cfg.ytick_count[j]),
+ )
+ elif not cfg.histogram[0]:
+ if min_v != max_v:
+ axis.set_ylim([min_v, max_v])
+ ylabels = np.linspace(min_v, max_v, int(cfg.ytick_count[j]))
+ if cfg.xlog[j] == "1":
+ axis.set_xscale("log")
+ else:
+ if tunit != "Dates":
+ if cfg.xformat[0]:
+ _set_formatted_ticks(axis, xlabels, cfg.xformat[j], "x")
+ elif not cfg.histogram[0]:
+ axis.set_xticks(xlabels)
+ if cfg.ylog[j] == "1":
+ axis.set_yscale("log")
+ else:
+ if cfg.yformat[0]:
+ _set_formatted_ticks(axis, ylabels, cfg.yformat[j], "y")
+ elif not cfg.histogram[0]:
+ axis.set_yticks(ylabels)
+ if cfg.legend_location[j] != "empty":
+ axis.legend(loc=cfg.legend_location[j])
+ if cfg.title[j] != "0" and cfg.hide_map_elements[3] == 0:
+ axis.set_title(cfg.title[j])
+ if cfg.remove_duplicate_labels and k + int(cfg.subplot_grid[1]) < len(
+ cfg.variables
+ ):
+ axis.tick_params(axis="x", which="both", bottom=False, labelbottom=False)
+ if len(cfg.variables) == len(cfg.cases[0]) and not cfg.subplot_grid[0]:
+ if cfg.csv:
+ generated_files.append(_save_summary_csv(cfg, deckn, var, quan, j))
+ return generated_files
+ generated_files.append(_save_summary_png(cfg, deckn, quan, j, fig))
+ return generated_files
+ if (
+ not cfg.subplot_grid[0] and len(cfg.variables) != len(cfg.cases[0])
+ ) or j == len(cfg.variables) - 1:
+ if (
+ len(cfg.legend_location) == j + 2
+ and j != 0
+ and len(axes.flat) - len(cfg.variables) > 0
+ ):
+ for jj, qua in enumerate(cfg.variables[: cfg.ncolors]):
+ for i, name in enumerate(cfg.cases[jj]):
+ time, var, tunit, vunit = read_series(
+ cfg,
+ name,
+ qua,
+ cfg.time_units[jj],
+ float(cfg.scale_factor[jj]),
+ i,
+ )
+ label = _get_label(cfg, name, jj, i)
+ if cfg.sensor or cfg.layer or cfg.distance[0]:
+ axes.flat[-1].plot(
+ time,
+ var,
+ color=cfg.colors[jj][i],
+ ls=cfg.linestyle[jj][i],
+ label=label,
+ lw=float(cfg.linewidth[jj][i]),
+ )
+ else:
+ axes.flat[-1].step(
+ time,
+ var,
+ color=cfg.colors[jj][i],
+ ls=cfg.linestyle[jj][i],
+ label=label,
+ lw=float(cfg.linewidth[jj][i]),
+ )
+ axes.flat[-1].axis("off")
+ axes.flat[-1].legend(loc=cfg.legend_location[-1])
+ for line in axes.flat[-1].get_lines():
+ line.remove()
+ for o in range(len(axes.flat) - len(cfg.variables) - 1):
+ fig.delaxes(axes.flat[-2 - o])
+ else:
+ for o in range(len(axes.flat) - len(cfg.variables)):
+ fig.delaxes(axes.flat[-1 - o])
+ generated_files.append(_save_summary_png(cfg, deckn, quan, j, fig))
+ plt.close()
+ return list(dict.fromkeys(generated_files))
+
+
+
+def _clean_name(name: str) -> str:
+ """Convert a variable expression to a filename-safe stem.
+
+ Parameters
+ ----------
+ name : str
+ Variable expression or proposed filename stem.
+
+ Returns
+ -------
+ str
+ Name with operators and separators replaced.
+
+ """
+ name = name.replace(" / ", "_over_")
+ name = name.replace(" ", "")
+ name = name.replace(":", "-")
+ return name
+
+
+def _get_deck_name(name: str) -> str:
+ """Get a display name from a case or include-file path.
+
+ Parameters
+ ----------
+ name : str
+ Case path or include filename.
+
+ Returns
+ -------
+ str
+ Lowercase basename without an ``.inc`` extension.
+
+ """
+ deckn = name.split("/")[-1].lower()
+ if ".inc" in deckn:
+ deckn = deckn[:-4]
+ return deckn
+
+
+def _get_label(cfg: PlopmConfig, name: str, var_index: int, name_index: int) -> str:
+ """Select the legend label for a plotted series.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized configuration containing ensemble cases and plot styles.
+ name : str
+ Simulation-case path.
+ var_index : int
+ Index of the plotted variable.
+ name_index : int
+ Index of the case within the variable group.
+
+ Returns
+ -------
+ str
+ User-defined label or a label derived from the case path.
+
+ """
+ label = name
+ if len(name.split("/")) > 1:
+ label = name.split("/")[-2] + "/" + name.split("/")[-1]
+ if cfg.legend_labels[0][0]:
+ label = cfg.legend_labels[var_index][name_index]
+ return label
+
+
+def _update_limits(
+ time: NDArray,
+ var: NDArray,
+ tunit: str,
+ min_t: float,
+ max_t: float,
+ min_v: float,
+ max_v: float,
+ xlow: float,
+ ylow: float,
+ first: bool,
+) -> tuple[float, float, float, float]:
+ """Update the data limits from one plotted series.
+
+ Parameters
+ ----------
+ time, var : np.ndarray
+ Time coordinates and variable values.
+ tunit : str
+ Time-axis label. ``"Dates"`` selects date handling.
+ min_t, max_t : float
+ Current time limits.
+ min_v, max_v : float
+ Current variable limits.
+ xlow, ylow : float
+ Lower bounds used to exclude invalid logarithmic values.
+ first : bool
+ Whether this is the first series included in the limits.
+
+ Returns
+ -------
+ tuple[float, float, float, float]
+ Updated ``(min_t, max_t, min_v, max_v)`` limits.
+
+ """
+ valid_var = var[var > ylow]
+ valid_time = time[time > xlow] if tunit != "Dates" else time
+ if valid_var.size == 0:
+ current_min_v = min_v if not first else 0
+ else:
+ current_min_v = np.min(valid_var)
+ current_max_v = np.nanmax(var) if np.any(~np.isnan(var)) else max_v
+ current_max_t = np.max(time)
+ if tunit != "Dates":
+ current_min_t = np.min(valid_time) if valid_time.size > 0 else min_t
+ else:
+ current_min_t = time[0]
+ if first:
+ return current_min_t, current_max_t, current_min_v, current_max_v
+ return (
+ min(min_t, current_min_t),
+ max(max_t, current_max_t),
+ min(min_v, current_min_v),
+ max(max_v, current_max_v),
+ )
+
+
+def _set_formatted_ticks(
+ axis: Axes, labels: NDArray, value_format: str, axis_name: str
+) -> None:
+ """Set explicitly formatted ticks on one axis.
+
+ Parameters
+ ----------
+ axis : matplotlib.axes.Axes
+ Axis to update.
+ labels : np.ndarray
+ Numeric tick locations.
+ value_format : str
+ Python format specification for each label.
+ axis_name : {"x", "y"}
+ Coordinate axis to update.
+
+ """
+ formatted_labels = [format(value, value_format) for value in labels]
+ ticks = [float(label) for label in formatted_labels]
+ if axis_name == "x":
+ axis.set_xticks(ticks)
+ axis.set_xticklabels(formatted_labels)
+ else:
+ axis.set_yticks(ticks)
+ axis.set_yticklabels(formatted_labels)
+
+
+def _save_summary_csv(
+ cfg: PlopmConfig, deckn: str, var: NDArray, quan: str, index: int
+) -> str:
+ """Write non-NaN summary values to a CSV file.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized configuration containing ensemble cases and plot styles.
+ deckn : str
+ Case name used in the default filename.
+ var : np.ndarray
+ Values to write.
+ quan : str
+ Variable expression used in the default filename.
+ index : int
+ Plot index used to select a custom filename.
+
+ Returns
+ -------
+ str
+ Name of the generated CSV file.
+
+ """
+ text = [f"{val}\n" for val in var if not np.isnan(val)]
+ name = _clean_name(f"{deckn}_{quan}")
+ if cfg.filename[index]:
+ name = cfg.filename[index]
+ filename = f"{name}.csv"
+ with open(
+ os.path.join(cfg.output_dir, filename),
+ "w",
+ encoding="utf8",
+ ) as file:
+ file.write("".join(text))
+ return filename
+
+
+def _save_summary_png(
+ cfg: PlopmConfig,
+ deckn: str,
+ quan: str,
+ index: int,
+ fig: Figure,
+) -> str:
+ """Save a summary figure as a PNG file.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized configuration containing ensemble cases and plot styles.
+ deckn : str
+ Case name used in the default filename.
+ quan : str
+ Variable expression used in the default filename.
+ index : int
+ Plot index used to select filename and resolution settings.
+ fig : matplotlib.figure.Figure
+ Figure to save.
+
+ Returns
+ -------
+ str
+ Name of the generated PNG file.
+
+ """
+ name = _clean_name(f"{deckn}_{quan}")
+ filename = f"{cfg.filename[index] if cfg.filename[index] else name}.png"
+ fig.savefig(
+ os.path.join(cfg.output_dir, filename),
+ bbox_inches="tight",
+ dpi=int(cfg.dpi[index]),
+ )
+ return filename
+
+
+def _plot_ensemble(
+ cfg: PlopmConfig, axes: Axes | np.ndarray
+) -> tuple[str, str, float, float, float, float]:
+ """Plot ensemble statistics for the first requested variable.
+
+ Each realization is interpolated to a shared coordinate array. Depending on
+ ``cfg.ensemble``, the function plots the mean, a one-standard-deviation band,
+ the bounding realizations, or both.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized configuration containing ensemble cases and plot styles.
+ axes : matplotlib.axes.Axes or np.ndarray
+ Axis, or array of axes, on which to draw the ensemble.
+
+ Returns
+ -------
+ tuple[str, str, float, float, float, float]
+ Time unit, value unit, and ``(min_t, max_t, min_v, max_v)`` limits.
+
+ """
+ axis = axes if isinstance(axes, Axes) else np.ravel(axes)[0]
+ thetime, timeeval = np.array([0]), np.array([0])
+ min_v, max_v = np.inf, -np.inf
+ hyst = 1
+ var_name = cfg.variables[0]
+ if (
+ var_name[:3] in ["krw", "krg"]
+ or var_name[:4]
+ in [
+ "krow",
+ "krog",
+ "pcow",
+ "pcog",
+ "pcwg",
+ ]
+ and var_name[-1] == "h"
+ ):
+ hyst = 2
+ for hyst_index in range(hyst):
+ for names_index, names in enumerate(cfg.cases):
+ label = cfg.case_labels[0][names_index] + " (mean)"
+ if len(label.split("/")) > 1:
+ label = label.split("/")[-2] + "/" + label.split("/")[-1]
+ if cfg.legend_labels[0][0]:
+ label = cfg.legend_labels[names_index][0]
+ tmp = []
+ for name_index, name in enumerate(names):
+ time, var, tunit, vunit = read_series(
+ cfg,
+ name,
+ var_name,
+ cfg.time_units[0],
+ float(cfg.scale_factor[0]),
+ name_index,
+ )
+ rng = int(1.0 * len(time) / hyst)
+ time = time[hyst_index * rng : (hyst_index + 1) * rng]
+ var = var[hyst_index * rng : (hyst_index + 1) * rng]
+ if time.size > thetime.size:
+ thetime = time.copy()
+ if tunit == "Dates":
+ time = np.array([value.timestamp() for value in time], dtype=float)
+ if time.size > timeeval.size:
+ timeeval = time.copy()
+ else:
+ timeeval = thetime
+ tmp.append(interp1d(time, var, bounds_error=False))
+ values = np.array([value(timeeval) for value in tmp])
+ with warnings.catch_warnings():
+ warnings.filterwarnings("ignore", message="Mean of empty slice")
+ warnings.filterwarnings("ignore", message="Degrees of freedom <= 0")
+ means = np.nanmean(values, axis=0)
+ stdev = np.nanstd(values, axis=0)
+ plot_label = label if hyst_index == hyst - 1 else None
+ axis.plot(
+ thetime,
+ means,
+ color=cfg.colors[0][names_index],
+ ls=cfg.linestyle[0][names_index],
+ label=plot_label,
+ lw=float(cfg.linewidth[0][names_index]),
+ )
+ if cfg.ensemble in [1, 3]:
+ if cfg.fill_between_style:
+ band_properties = cfg.fill_between_style.split(",")
+ color = band_properties[2 * names_index]
+ alpha = float(band_properties[2 * names_index + 1])
+ else:
+ color = cfg.colors[0][names_index]
+ alpha = 0.2
+ lower_band = means - stdev
+ upper_band = means + stdev
+ axis.fill_between(
+ thetime, lower_band, upper_band, color=color, alpha=alpha
+ )
+ if np.any(~np.isnan(lower_band)):
+ min_v = min(min_v, np.nanmin(lower_band))
+ if np.any(~np.isnan(upper_band)):
+ max_v = max(max_v, np.nanmax(upper_band))
+ if cfg.ensemble in [2, 3]:
+ ensemble_index = len(cfg.cases) + names_index
+ maxs = np.nansum(values + means, axis=1)
+ mins = np.nansum(values - means, axis=1)
+ maxs = np.where(maxs == np.max(maxs))[0][0]
+ mins = np.where(mins == np.min(mins))[0][0]
+ labell = names[mins] + " (lower)"
+ labelu = names[maxs] + " (upper)"
+ if cfg.legend_labels[0][0]:
+ labell = cfg.legend_labels[names_index][1]
+ labelu = cfg.legend_labels[names_index][2]
+ lower_label = labell if hyst_index == hyst - 1 else None
+ upper_label = labelu if hyst_index == hyst - 1 else None
+ axis.plot(
+ thetime,
+ values[mins],
+ color=cfg.colors[0][ensemble_index],
+ ls=cfg.linestyle[0][ensemble_index],
+ label=lower_label,
+ lw=float(cfg.linewidth[0][names_index]),
+ )
+ axis.plot(
+ thetime,
+ values[maxs],
+ color=cfg.colors[0][ensemble_index],
+ ls=cfg.linestyle[0][ensemble_index],
+ label=upper_label,
+ lw=float(cfg.linewidth[0][names_index]),
+ )
+ if np.any(~np.isnan(values[mins])):
+ min_v = min(min_v, np.nanmin(values[mins]))
+ if np.any(~np.isnan(values[maxs])):
+ max_v = max(max_v, np.nanmax(values[maxs]))
+ min_t, max_t = thetime[0], thetime[-1]
+ return tunit, vunit, min_t, max_t, min_v, max_v
+
+
+
+
+
+# SPDX-FileCopyrightText: 2024-2026 NORCE Research AS
+# SPDX-License-Identifier: GPL-3.0
+# pylint: disable=W3301,W0123,R0912,R0915,R0914,R1702,W0611,R0913,R0917,C0302,C0115,R0916,E1102
+
+"""Create two-dimensional maps and animations from OPM results.
+
+The module prepares grid geometry, maps three-dimensional properties onto
+selected slices, and writes PNG or GIF output with optional masks, differences,
+well and fault overlays, and shared color limits.
+"""
+
+import datetime
+import sys
+from collections.abc import Iterable
+from contextlib import nullcontext
+from typing import Any
+
+import colorcet # noqa: F401 # registers colorcet colormaps with matplotlib
+import matplotlib
+import matplotlib.pyplot as plt
+import matplotlib.ticker as mticker
+import numpy as np
+from alive_progress import alive_bar
+from matplotlib import animation, colors
+from matplotlib.animation import FuncAnimation, writers
+from matplotlib.artist import Artist
+from matplotlib.axes import Axes
+from matplotlib.cm import ScalarMappable
+from matplotlib.figure import Figure
+from matplotlib.ticker import LogFormatter
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+from mpl_toolkits.axes_grid1.axes_divider import AxesDivider
+from numpy.typing import NDArray
+
+from plopm.config.config import PlopmConfig, SimData
+from plopm.utils.mapping import (
+ get_xy_slice,
+ get_xz_slice,
+ get_yz_slice,
+ map_xy,
+ map_xz,
+ map_yz,
+ transform_grid,
+)
+from plopm.utils.readers import (
+ get_faults,
+ get_wells,
+ read_case,
+ read_csv_grid,
+ read_quantity,
+ time_unit,
+)
+from plopm.utils.terminal import cli_error_value, plopm_error
+
+
+
+[docs]
+def make_maps(cfg: PlopmConfig) -> list[str]:
+ """Create the requested spatial maps and animations.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map configuration.
+
+ Returns
+ -------
+ list[str]
+ Names of the generated PNG and GIF files.
+
+ """
+ generated_files: list[str] = []
+ skip = 0
+ if (
+ cfg.subplot_grid[0]
+ and len(cfg.variables) > 1
+ and len(cfg.restart) == 1
+ and len(cfg.cases[0]) == 1
+ ):
+ skip = 1
+ if cfg.subplot_grid[0]:
+ fig, axis = _create_figure(int(cfg.subplot_grid[0]), int(cfg.subplot_grid[1]))
+ sub1 = int(cfg.subplot_grid[1])
+ else:
+ fig, axis = _create_figure(1, 1, "compressed")
+ sub1 = 1
+ if cfg.subplot_grid[0] and cfg.gif and len(cfg.cases[0]) > 1:
+ _, _, _, cmin, cmax, diffa = _get_clim(cfg)
+ maska = _get_masks(cfg) if cfg.mask_variable else []
+ deckd = _case_name(cfg.difference_input) if cfg.difference_input else ""
+ fig, axis = _create_figure(
+ int(cfg.subplot_grid[0]), int(cfg.subplot_grid[1]), "compressed"
+ )
+ axes = _normalize_axis(axis)
+ data, xc, yc, named, slice_title, slice_name, mx, my, xname, yname = (
+ _prepare_map(cfg, cfg.cases[0][0], 0)
+ )
+ original_loc, cb = _prepare_colorbars(axes)
+ _delete_extra_axes(axes, len(cfg.cases[0]), fig)
+ im_ani = animation.FuncAnimation(
+ fig,
+ _draw_frame,
+ fargs=(
+ cfg.cases[0][0],
+ fig,
+ axes,
+ original_loc,
+ cb,
+ cmin,
+ cmax,
+ maska,
+ diffa,
+ named,
+ deckd,
+ slice_title,
+ slice_name,
+ cfg,
+ generated_files,
+ 0,
+ data,
+ xc,
+ yc,
+ skip,
+ sub1,
+ mx,
+ my,
+ xname,
+ yname,
+ ),
+ frames=len(data.steps),
+ interval=cfg.gif_interval,
+ blit=False,
+ repeat=False,
+ )
+ generated_files.append(
+ _save_animation(
+ cfg, im_ani, cfg.filename[0] if cfg.filename[0] else cfg.variables[0]
+ )
+ )
+ elif cfg.subplot_grid[0] and cfg.gif and len(cfg.variables) > 1:
+ data, xc, yc, cmin, cmax, diffa = _get_clim(cfg)
+ deckd = _case_name(cfg.difference_input) if cfg.difference_input else ""
+ data, xc, yc, named, slice_title, slice_name, mx, my, xname, yname = (
+ _prepare_map(cfg, cfg.cases[0][0], 0)
+ )
+ maska = _get_masks(cfg) if cfg.mask_variable else []
+ if len(data.steps) > 1:
+ fig, axis = _create_figure(
+ int(cfg.subplot_grid[0]), int(cfg.subplot_grid[1])
+ )
+ axes = _normalize_axis(axis)
+ plt.tight_layout(pad=1.7)
+ original_loc, cb = _prepare_colorbars(axes)
+ _delete_extra_axes(axes, len(cfg.variables), fig)
+ im_ani = animation.FuncAnimation(
+ fig,
+ _draw_frame,
+ fargs=(
+ cfg.cases[0][0],
+ fig,
+ axes,
+ original_loc,
+ cb,
+ cmin,
+ cmax,
+ maska,
+ diffa,
+ named,
+ deckd,
+ slice_title,
+ slice_name,
+ cfg,
+ generated_files,
+ 0,
+ data,
+ xc,
+ yc,
+ skip,
+ sub1,
+ mx,
+ my,
+ xname,
+ yname,
+ ),
+ frames=len(data.steps),
+ interval=cfg.gif_interval,
+ blit=False,
+ repeat=False,
+ )
+ generated_files.append(
+ _save_animation(cfg, im_ani, cfg.filename[0] if cfg.filename[0] else named)
+ )
+ else:
+ _, _, _, cmin, cmax, diffa = _get_clim(cfg)
+ maska = _get_masks(cfg) if cfg.mask_variable else []
+ deckd = _case_name(cfg.difference_input) if cfg.difference_input else ""
+ data, xc, yc, named, slice_title, slice_name, mx, my, xname, yname = (
+ _prepare_map(cfg, cfg.cases[0][0], 0)
+ )
+ for n, var in enumerate(cfg.variables):
+ if len(data.steps) > 1:
+ if cfg.subplot_grid[0]:
+ fig, axis = _create_figure(
+ int(cfg.subplot_grid[0]), int(cfg.subplot_grid[1])
+ )
+ else:
+ fig, axis = _create_figure(1, 1)
+ if not cfg.subplot_grid[0] and not cfg.gif:
+ plt.close()
+ fig, axis = _create_figure(1, 1, "tight")
+ axes = _normalize_axis(axis)
+ original_loc, cb = _prepare_colorbars(axes)
+ if len(data.steps) > 1:
+ _delete_extra_axes(axes, len(data.steps), fig)
+ if cfg.gif and len(data.steps) > 1:
+ im_ani = animation.FuncAnimation(
+ fig,
+ _draw_frame,
+ fargs=(
+ cfg.cases[0][0],
+ fig,
+ axes,
+ original_loc,
+ cb,
+ cmin,
+ cmax,
+ maska,
+ diffa,
+ named,
+ deckd,
+ slice_title,
+ slice_name,
+ cfg,
+ generated_files,
+ n,
+ data,
+ xc,
+ yc,
+ skip,
+ sub1,
+ mx,
+ my,
+ xname,
+ yname,
+ ),
+ frames=len(data.steps),
+ interval=cfg.gif_interval,
+ blit=False,
+ repeat=False,
+ )
+ name = f"{cfg.filename[0] if cfg.filename[0] else named + '_' + var}"
+ generated_files.append(_save_animation(cfg, im_ani, name))
+ else:
+ if len(cfg.cases[0]) > 1:
+ _delete_extra_axes(axes, len(cfg.cases[0]), fig)
+ if len(data.steps) > 1 and len(cfg.cases[0]) == len(data.steps):
+ if not cfg.subplot_grid[0]:
+ fig, axis = _create_figure(1, 1)
+ axes = _normalize_axis(axis)
+ original_loc, cb = _prepare_colorbars(axes)
+ _draw_frame(
+ 0,
+ cfg.cases[0][0],
+ fig,
+ axes,
+ original_loc,
+ cb,
+ cmin,
+ cmax,
+ maska,
+ diffa,
+ named,
+ deckd,
+ slice_title,
+ slice_name,
+ cfg,
+ generated_files,
+ n,
+ data,
+ xc,
+ yc,
+ skip,
+ sub1,
+ mx,
+ my,
+ xname,
+ yname,
+ )
+ else:
+ for t, _ in enumerate(data.steps):
+ if not cfg.subplot_grid[0]:
+ plt.close()
+ fig, axis = _create_figure(1, 1)
+ axes = _normalize_axis(axis)
+ original_loc, cb = _prepare_colorbars(axes)
+ _draw_frame(
+ t,
+ cfg.cases[0][0],
+ fig,
+ axes,
+ original_loc,
+ cb,
+ cmin,
+ cmax,
+ maska,
+ diffa,
+ named,
+ deckd,
+ slice_title,
+ slice_name,
+ cfg,
+ generated_files,
+ n,
+ data,
+ xc,
+ yc,
+ skip,
+ sub1,
+ mx,
+ my,
+ xname,
+ yname,
+ )
+ return list(dict.fromkeys(generated_files))
+
+
+
+def _prepare_map(
+ cfg: PlopmConfig, deck: str, n: int
+) -> tuple[SimData, NDArray, NDArray, str, str, str, int, int, str, str]:
+ """Prepare simulation data and coordinates for one map.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map configuration.
+ deck : str
+ Simulation-case stem or CSV input path.
+ n : int
+ Case or map index used to select configuration values.
+
+ Returns
+ -------
+ tuple
+ Simulation data, coordinate meshes, case and slice labels, mesh
+ dimensions, and coordinate-axis names.
+
+ """
+ if cfg.csv_columns[n][0]:
+ xc, yc, mx, my, xname, yname = read_csv_grid(cfg, deck, n)
+ slice_title, slice_name = "", ""
+ data = SimData(steps=cfg.restart)
+ else:
+ data = read_case(
+ deck, cfg.gif, cfg.vtk, cfg.variables, cfg.restart, cfg.filters, n
+ )
+ slice_value = cfg.slice[n]
+ if slice_value[0][0] != -2:
+ xc, yc, slice_title, slice_name, mx, my, xname, yname = get_yz_slice(
+ cfg, data, n
+ )
+ elif slice_value[1][0] != -2:
+ xc, yc, slice_title, slice_name, mx, my, xname, yname = get_xz_slice(
+ cfg, data, n
+ )
+ else:
+ xc, yc, slice_title, slice_name, mx, my, xname, yname = get_xy_slice(
+ cfg, data, n
+ )
+ if int(cfg.rotation[n]) != 0 or cfg.translation[n] != ["[0", "0]"]:
+ xc, yc = transform_grid(cfg, n, xc, yc)
+ return (
+ data,
+ xc,
+ yc,
+ deck.rsplit("/", 1)[-1].lower(),
+ slice_title,
+ slice_name,
+ mx,
+ my,
+ xname,
+ yname,
+ )
+
+
+def _create_figure(
+ rows: int = 1,
+ columns: int = 1,
+ layout: str | None = None,
+) -> tuple[Figure, Axes]:
+ """Create a Matplotlib figure and axes.
+
+ Parameters
+ ----------
+ rows, columns : int, default: 1
+ Number of subplot rows and columns.
+ layout : str, optional
+ Matplotlib layout engine.
+
+ Returns
+ -------
+ tuple
+ Created figure and axes.
+
+ """
+ plt.close()
+ if layout:
+ fig, axes = plt.subplots(rows, columns, layout=layout)
+ else:
+ fig, axes = plt.subplots(rows, columns)
+ return fig, axes
+
+
+def _normalize_axis(axes: Axes | NDArray[Any]) -> NDArray:
+ """Return axes as a one-dimensional-compatible array.
+
+ Parameters
+ ----------
+ axes : matplotlib.axes.Axes or np.ndarray
+ Axes returned by Matplotlib.
+
+ Returns
+ -------
+ np.ndarray
+ Array containing the supplied axes.
+
+ """
+ if isinstance(axes, np.ndarray):
+ return axes
+ return np.array([axes])
+
+
+def _prepare_colorbars(axes: NDArray[Any]) -> tuple[list[Any], list[str]]:
+ """Initialize colorbar state for each subplot.
+
+ Parameters
+ ----------
+ axes : np.ndarray
+ Subplot axes.
+
+ Returns
+ -------
+ tuple[list, list]
+ Original axes locators and empty colorbar slots.
+
+ """
+ original_loc, cb = [], []
+ for axis in axes.flat:
+ original_loc.append(axis.get_axes_locator())
+ cb.append("")
+ return original_loc, cb
+
+
+def _delete_extra_axes(axes: NDArray[Any], keep: int, fig: Figure) -> None:
+ """Remove unused subplot axes.
+
+ Parameters
+ ----------
+ axes : np.ndarray
+ Subplot axes.
+ keep : int
+ Number of axes to retain.
+ fig : matplotlib.figure.Figure
+ Figure containing the axes.
+
+ """
+ for o in range(max(0, len(axes.flat) - keep)):
+ axis_to_remove = axes.flat[-1 - o]
+ if axis_to_remove in fig.axes:
+ fig.delaxes(axis_to_remove)
+
+
+def _save_animation(cfg: PlopmConfig, im_ani: FuncAnimation, name: str) -> str:
+ """Save a Matplotlib animation as a GIF.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map configuration.
+ im_ani : matplotlib.animation.FuncAnimation
+ Animation to save.
+ name : str
+ Output filename without extension.
+
+ Returns
+ -------
+ str
+ Name of the generated GIF file.
+
+ """
+ filename = f"{name}.gif"
+ output_path = f"{cfg.output_dir}/{filename}"
+ if cfg.gif_loop or not writers.is_available("ffmpeg"):
+ im_ani.save(output_path)
+ else:
+ im_ani.save(output_path, extra_args=["-loop", "-1"])
+ return filename
+
+
+def _map_values(
+ cfg: PlopmConfig,
+ data: SimData,
+ var: str,
+ values: NDArray,
+ slice_index: int,
+ map_index: int,
+ mx: int,
+ my: int,
+ use_csv: bool = False,
+) -> NDArray:
+ """Map quantity values onto the selected two-dimensional slice.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map configuration.
+ data : SimData
+ Loaded simulation data.
+ var : str
+ Variable name.
+ values : np.ndarray
+ Values in active-cell or CSV order.
+ slice_index, map_index : int
+ Indices selecting the slice and its mapping settings.
+ mx, my : int
+ Mapped grid dimensions.
+ use_csv : bool, default: False
+ Whether values already use the two-dimensional CSV layout.
+
+ Returns
+ -------
+ np.ndarray
+ Values arranged on the selected map.
+
+ """
+ if use_csv:
+ quaa = np.asarray(values).copy()
+ elif cfg.slice[slice_index][0][0] != -2:
+ quaa = map_yz(cfg, data, var, values, map_index, mx, my)
+ elif cfg.slice[slice_index][1][0] != -2:
+ quaa = map_xz(cfg, data, var, values, map_index, mx, my)
+ else:
+ quaa = map_xy(cfg, data, var, values, map_index, mx, my)
+ return quaa
+
+
+def _get_clim(
+ cfg: PlopmConfig,
+) -> tuple[SimData, NDArray, NDArray, list[float], list[float], list[NDArray]]:
+ """Determine color limits and cached difference maps.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map configuration.
+
+ Returns
+ -------
+ tuple
+ Last loaded simulation data, coordinate meshes, color minima and
+ maxima, and cached difference arrays.
+
+ """
+ cmin, cmax = [float("inf")], [float("-inf")]
+ diffa: list[NDArray] = []
+ xc, yc = np.empty(0), np.empty(0)
+ if (cfg.rst_range and cfg.png and not cfg.subplot_grid[0]) or (
+ cfg.clim[0][0] and not cfg.difference_input
+ ):
+ return SimData(), xc, yc, cmin, cmax, diffa
+
+ if cfg.restart[0] == -1 and cfg.gif:
+ data = read_case(
+ cfg.cases[0][0], cfg.gif, cfg.vtk, cfg.variables, cfg.restart, cfg.filters
+ )
+ else:
+ data = SimData(steps=cfg.restart)
+ if cfg.difference_input:
+ var = cfg.variables[0]
+ for t, _ in enumerate(data.steps):
+ data, xc, yc, _, _, _, mx, my, _, _ = _prepare_map(
+ cfg, cfg.difference_input, 1
+ )
+ _, values = read_quantity(
+ cfg.difference_input,
+ data,
+ var,
+ data.steps[t],
+ float(cfg.scale_factor[0]),
+ cfg.mass_vars,
+ cfg.mass_vars + cfg.mass_fracs,
+ cfg.caprock_vars,
+ cfg.stress_coefficient,
+ cfg.filters[0],
+ cfg.gif,
+ cfg.min_threshold[0],
+ cfg.max_threshold[0],
+ cfg.csv_columns[0],
+ )
+ quaa = _map_values(cfg, data, var, values, 1, 1, mx, my)
+ diffa.append(quaa.copy())
+ if len(cfg.variables) == len(cfg.cases[0]) and len(cfg.cases[0]) > 1:
+ for m, var in enumerate(cfg.variables):
+ cmin.append(cmin[-1])
+ cmax.append(cmax[-1])
+ for t, _ in enumerate(data.steps):
+ data, xc, yc, _, _, _, mx, my, _, _ = _prepare_map(
+ cfg, cfg.cases[0][m], m
+ )
+ _, values = read_quantity(
+ cfg.cases[0][m],
+ data,
+ var,
+ data.steps[t],
+ float(cfg.scale_factor[m]),
+ cfg.mass_vars,
+ cfg.mass_vars + cfg.mass_fracs,
+ cfg.caprock_vars,
+ cfg.stress_coefficient,
+ cfg.filters[0],
+ cfg.gif,
+ cfg.min_threshold[m],
+ cfg.max_threshold[m],
+ cfg.csv_columns[0],
+ )
+ quaa = _map_values(cfg, data, var, values, m, m, mx, my)
+ _apply_diff_and_log(cfg, diffa, quaa, m, t)
+ _update_color_range(quaa, cmin, cmax)
+ else:
+ for m, var in enumerate(cfg.variables):
+ cmin.append(cmin[-1])
+ cmax.append(cmax[-1])
+ for n, deck in enumerate(cfg.cases[0]):
+ for t, _ in enumerate(data.steps):
+ data, xc, yc, _, _, _, mx, my, _, _ = _prepare_map(cfg, deck, n)
+ _, values = read_quantity(
+ deck,
+ data,
+ var,
+ data.steps[t],
+ float(cfg.scale_factor[m]),
+ cfg.mass_vars,
+ cfg.mass_vars + cfg.mass_fracs,
+ cfg.caprock_vars,
+ cfg.stress_coefficient,
+ cfg.filters[n],
+ cfg.gif,
+ cfg.min_threshold[m],
+ cfg.max_threshold[m],
+ cfg.csv_columns[n],
+ )
+ quaa = _map_values(
+ cfg, data, var, values, n, n, mx, my, cfg.csv_columns[n][0]
+ )
+ _apply_diff_and_log(cfg, diffa, quaa, m, t)
+ _update_color_range(quaa, cmin, cmax)
+ return data, xc, yc, cmin, cmax, diffa
+
+
+def _apply_diff_and_log(
+ cfg: PlopmConfig,
+ diffa: list[NDArray],
+ quaa: NDArray,
+ var_index: int,
+ restart_index: int,
+) -> None:
+ """Apply difference and logarithmic transformations in place.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map configuration.
+ diffa : list[np.ndarray]
+ Mapped difference values.
+ quaa : np.ndarray
+ Mapped values to transform.
+ var_index : int
+ Variable index.
+ restart_index : int
+ Restart-step index used to select a cached difference map.
+
+ """
+ if cfg.difference_input:
+ quaa -= diffa[restart_index]
+ if int(cfg.color_log[var_index]) == 1:
+ quaa[quaa <= 0] = np.nan
+
+
+def _update_color_range(quaa: NDArray, cmin: list[float], cmax: list[float]) -> None:
+ """Update the current finite color range.
+
+ Parameters
+ ----------
+ quaa : np.ndarray
+ Mapped values included in the color range.
+ cmin : list[float]
+ Color minima.
+ cmax : list[float]
+ Color minima.
+
+ """
+ if np.any(~np.isnan(quaa)):
+ cmin[-2] = min(cmin[-2], np.nanmin(quaa))
+ cmax[-2] = max(cmax[-2], np.nanmax(quaa))
+
+
+def _get_masks(cfg: PlopmConfig) -> list[NDArray]:
+ """Read and map masks for all configured cases.
+
+ Parameters
+ ----------
+ cfg : PlopmConfig
+ Initialized map configuration.
+
+ Returns
+ -------
+ list[np.ndarray]
+ Mapped mask arrays.
+
+ """
+ maska = []
+ var = cfg.mask_variable
+ for n, deck in enumerate(cfg.cases[0]):
+ data, _, _, _, _, _, mx, my, _, _ = _prepare_map(cfg, deck, n)
+ _, values = read_quantity(
+ deck,
+ data,
+ var,
+ 0,
+ float(cfg.scale_factor[0]),
+ cfg.mass_vars,
+ cfg.mass_vars + cfg.mass_fracs,
+ cfg.caprock_vars,
+ cfg.stress_coefficient,
+ cfg.filters[n],
+ cfg.gif,
+ cfg.min_threshold[0],
+ cfg.max_threshold[0],
+ cfg.csv_columns[n],
+ )
+ maska.append(_map_values(cfg, data, var, values, n, n, mx, my))
+ return maska
+
+
+def _case_name(deck: str) -> str:
+ """Get a lowercase case name from a path.
+
+ Parameters
+ ----------
+ deck : str
+ Simulation-case path.
+
+ Returns
+ -------
+ str
+ Final path component in lowercase.
+
+ """
+ if len(deck.split("/")) > 1:
+ return deck.split("/")[-1].lower()
+ return deck.lower()
+
+
+def _draw_frame(
+ t: int,
+ deck: str,
+ fig: Figure,
+ axes: Any,
+ original_loc: list[Any],
+ cb: list[str],
+ cmin: list[float],
+ cmax: list[float],
+ maska: list[Any],
+ diffa: list[NDArray],
+ named: str,
+ deckd: str,
+ slice_title: str,
+ slice_name: str,
+ cfg: PlopmConfig,
+ generated_files: list[str],
+ n: int,
+ data: SimData,
+ xc: NDArray,
+ yc: NDArray,
+ skip: int,
+ sub1: int,
+ mx: int,
+ my: int,
+ xname: str,
+ yname: str,
+) -> Iterable[Artist]:
+ """Draw all maps belonging to one animation frame.
+
+ This dispatcher selects cases, variables, restart steps, and subplot
+ positions before delegating each map to :func:`draw_map`.
+
+ Parameters
+ ----------
+ t : int
+ Animation-frame or restart-step index.
+ deck : str
+ Primary simulation-case stem.
+ fig : matplotlib.figure.Figure
+ Figure receiving the maps.
+ axes : matplotlib.axes.Axes or np.ndarray
+ Target axes.
+ original_loc, cb : list
+ Original axes locators and active colorbars.
+ cmin, cmax : list[float]
+ Color limits for each variable or map.
+ maska, diffa : list
+ Mapped masks and cached difference arrays.
+ named, deckd : str
+ Display names for the primary and difference cases.
+ slice_title, slice_name : str
+ Human-readable slice descriptions.
+ cfg : PlopmConfig
+ Initialized map configuration.
+ generated_files : list[str]
+ Generated filenames updated during rendering.
+ n : int
+ Current variable or case index.
+ data : SimData
+ Loaded simulation data.
+ xc, yc : np.ndarray
+ Coordinate meshes.
+ skip, sub1 : int
+ Subplot-control values.
+ mx, my : int
+ Mapped grid dimensions.
+ xname, yname : str
+ Coordinate-axis names.
+
+ Returns
+ -------
+ list[matplotlib.artist.Artist]
+ Empty artist list required by the animation callback.
+
+ """
+ k = t
+ if not cfg.subplot_grid[0]:
+ k = 0
+ elif len(data.steps) == 1:
+ k = n
+ if cfg.subplot_grid[0] and len(cfg.cases[0]) > 1:
+ show_progress = sys.stdout.isatty()
+ if show_progress:
+ bar_ctx = alive_bar(len(cfg.cases[0]), bar="fish")
+ else:
+ bar_ctx = nullcontext()
+ with bar_ctx as bar_animation:
+ if len(cfg.variables) > 1:
+ cmax = [np.max(cmax)] * len(cmax)
+ cmin = [np.min(cmin)] * len(cmin)
+ for nn, deckl in enumerate(cfg.cases[0]):
+ if show_progress:
+ bar_animation()
+ (
+ data,
+ xc,
+ yc,
+ named,
+ slice_title,
+ slice_name,
+ mx,
+ my,
+ xname,
+ yname,
+ ) = _prepare_map(cfg, deckl, nn)
+ _draw_map(
+ deckl,
+ fig,
+ axes,
+ original_loc,
+ cb,
+ cmin,
+ cmax,
+ maska,
+ diffa,
+ named,
+ deckd,
+ slice_title,
+ slice_name,
+ cfg,
+ generated_files,
+ data,
+ t,
+ nn,
+ nn,
+ xc,
+ yc,
+ sub1,
+ mx,
+ my,
+ xname,
+ yname,
+ )
+ else:
+ for nn, deckl in enumerate(cfg.cases[0]):
+ if show_progress:
+ bar_animation()
+ (
+ data,
+ xc,
+ yc,
+ named,
+ slice_title,
+ slice_name,
+ mx,
+ my,
+ xname,
+ yname,
+ ) = _prepare_map(cfg, deckl, nn)
+ if len(data.steps) > 1 and len(cfg.cases[0]) == len(data.steps):
+ _draw_map(
+ deckl,
+ fig,
+ axes,
+ original_loc,
+ cb,
+ cmin,
+ cmax,
+ maska,
+ diffa,
+ named,
+ deckd,
+ slice_title,
+ slice_name,
+ cfg,
+ generated_files,
+ data,
+ nn,
+ 0,
+ nn,
+ xc,
+ yc,
+ sub1,
+ mx,
+ my,
+ xname,
+ yname,
+ )
+ else:
+ _draw_map(
+ deckl,
+ fig,
+ axes,
+ original_loc,
+ cb,
+ cmin,
+ cmax,
+ maska,
+ diffa,
+ named,
+ deckd,
+ slice_title,
+ slice_name,
+ cfg,
+ generated_files,
+ data,
+ t,
+ 0,
+ nn,
+ xc,
+ yc,
+ sub1,
+ mx,
+ my,
+ xname,
+ yname,
+ )
+ elif cfg.subplot_grid[0] and len(cfg.variables) > 1 and skip == 0:
+ show_progress = sys.stdout.isatty()
+ if show_progress:
+ bar_ctx = alive_bar(len(cfg.variables), bar="fish")
+ else:
+ bar_ctx = nullcontext()
+ with bar_ctx as bar_animation:
+ for nn, _ in enumerate(cfg.variables):
+ if show_progress:
+ bar_animation()
+ _draw_map(
+ deck,
+ fig,
+ axes,
+ original_loc,
+ cb,
+ cmin,
+ cmax,
+ maska,
+ diffa,
+ named,
+ deckd,
+ slice_title,
+ slice_name,
+ cfg,
+ generated_files,
+ data,
+ t,
+ nn,
+ nn,
+ xc,
+ yc,
+ sub1,
+ mx,
+ my,
+ xname,
+ yname,
+ )
+ else:
+ _draw_map(
+ deck,
+ fig,
+ axes,
+ original_loc,
+ cb,
+ cmin,
+ cmax,
+ maska,
+ diffa,
+ named,
+ deckd,
+ slice_title,
+ slice_name,
+ cfg,
+ generated_files,
+ data,
+ t,
+ n,
+ k,
+ xc,
+ yc,
+ sub1,
+ mx,
+ my,
+ xname,
+ yname,
+ )
+ return []
+
+
+def _set_axis(
+ fig: Figure,
+ axes: Any,
+ cfg: PlopmConfig,
+ data: SimData,
+ name: str,
+ n: int,
+ t: int,
+ k: int,
+ n_s: int,
+ unit: str,
+ xc: NDArray,
+ yc: NDArray,
+ extinf: float,
+ named: str,
+ deckd: str,
+ defcol: bool,
+ slice_title: str,
+ feature_id: int,
+) -> None:
+ """Configure labels, limits, ticks, and annotations for a map axis.
+
+ Parameters
+ ----------
+ fig : matplotlib.figure.Figure
+ Figure containing the map.
+ axes : matplotlib.axes.Axes or np.ndarray
+ Map axes.
+ cfg : PlopmConfig
+ Initialized map configuration.
+ data : SimData
+ Loaded simulation data.
+ name : str
+ Variable name.
+ n, t, k, n_s : int
+ Variable, restart, subplot, and slice indices.
+ unit : str
+ Variable unit label.
+ xc, yc : np.ndarray
+ Coordinate meshes.
+ extinf : float
+ Padding added to map extents.
+ named, deckd : str
+ Display names for the primary and difference cases.
+ defcol : bool
+ Whether default categorical colors are used.
+ slice_title : str
+ Human-readable slice description.
+ feature_id : int
+ Number assigned to the active well or fault feature.
+
+ """
+ unrst = data.unrst
+ nx = data.nx
+ ny = data.ny
+ nz = data.nz
+ actind = data.active_idx
+ restart = data.steps
+ porv = data.active_pv
+ axis = axes.flat[k]
+ name_lower = name.lower()
+ is_discrete_num = (
+ "num" in name
+ and (cfg.colormaps[n] in cfg.disc_colormaps or defcol)
+ and cfg.discrete
+ )
+ namet, time = name, ""
+ if cfg.time_units[0] == "dates":
+ date_values = unrst["INTEHEAD", restart[t]]
+ date = datetime.date(
+ date_values[66],
+ date_values[65],
+ date_values[64],
+ )
+ time = f" {date}"
+ elif cfg.time_units[0] == "empty":
+ pass
+ else:
+ tskl, tunit = time_unit(cfg.time_units[0])
+ tunit = tunit[5:]
+ if unrst and unrst.count("DOUBHEAD", 0):
+ time = f" {tskl*unrst['DOUBHEAD', restart[t]][0]:.0f} {tunit}"
+ elif cfg.time_units[0] in ["s", "m", "h", "d", "w", "y"]:
+ time = f" {restart[t]:.0f} {tunit}"
+ else:
+ time = f" {restart[t]:.0f} [{cfg.time_units[0]}]"
+ if cfg.equal_aspect:
+ axis.axis("scaled")
+ extra = ""
+ if name_lower == "porv":
+ extra = f", sum={np.sum(porv):.3e}"
+ elif name_lower in cfg.mass_vars and cfg.difference_input:
+ extra = f", |sum|={extinf:.3e} {unit}"
+ elif name_lower in cfg.mass_vars:
+ extra = f", sum={extinf:.3e} {unit}"
+ elif cfg.difference_input:
+ extra = f", |sum|={extinf:.3e}"
+ elif cfg.variables[0] in ["wells", "faults"]:
+ time = ""
+ namet = f"Total no. {name} = {feature_id-1}, "
+ elif is_discrete_num:
+ time = ""
+ namet = ""
+ if cfg.csv_columns[n][0]:
+ tslice = ""
+ elif cfg.variables[0] in ["wells", "faults"] or is_discrete_num:
+ tslice = slice_title[2:]
+ else:
+ tslice = slice_title
+ if (
+ cfg.subplot_grid[0]
+ and len(cfg.cases[0]) > 1
+ and cfg.title[k] == "0"
+ and cfg.hide_map_elements[3] == 0
+ ):
+ if name_lower == "porv":
+ named += f" (total porv={np.sum(data.porv)})"
+ axis.set_title(named)
+ if k == 0 and cfg.suptitle != "0":
+ fig.suptitle(f"{time[1:]}")
+ elif cfg.subplot_grid[0] and len(cfg.variables) > 1 and cfg.title[k] == "0":
+ if k == 0 and cfg.suptitle != "0":
+ fig.suptitle(f"{named}{time}")
+ elif (
+ cfg.gif
+ and len(cfg.variables) == 1
+ and cfg.title[k] == "0"
+ and cfg.hide_map_elements[3] == 0
+ ):
+ if cfg.difference_input:
+ axis.set_title(f"{named}-{deckd}{time}")
+ else:
+ axis.set_title(f"{named}{time}")
+ elif (
+ cfg.gif
+ and len(cfg.variables) == 1
+ and cfg.title[k] != "0"
+ and cfg.hide_map_elements[3] == 0
+ ):
+ if not cfg.csv_columns[n][0]:
+ axis.set_title(f"{cfg.title[k]} {time}")
+ else:
+ axis.set_title(f"{cfg.title[k]}")
+ fig.suptitle(time)
+ elif (
+ len(restart) > 1
+ and cfg.subplot_grid[0]
+ and len(cfg.cases[0]) == 1
+ and cfg.title[k] == "0"
+ and cfg.hide_map_elements[3] == 0
+ ):
+ axis.set_title(f"{unrst['DOUBHEAD', restart[t]][0]} days")
+ if k == 0 and cfg.suptitle != "0":
+ if cfg.difference_input:
+ fig.suptitle(f"{named}-{deckd}")
+ else:
+ fig.suptitle(f"{named}")
+ elif cfg.hide_map_elements[3] == 0 and cfg.title[k] == "0":
+ if cfg.difference_input:
+ axis.set_title(f"{named}-{deckd}" + tslice + extra + time)
+ else:
+ axis.set_title(namet + tslice + extra + time)
+ elif cfg.subplot_grid[0] and len(cfg.cases[0]) > 1:
+ if k == 0 and cfg.suptitle != "0":
+ if cfg.gif and cfg.csv_columns[n][0]:
+ fig.suptitle(f"{restart[t]} {cfg.time_units[0]}")
+ elif unrst:
+ fig.suptitle(f"{unrst['DOUBHEAD', restart[t]][0]} days")
+ else:
+ fig.suptitle(f"{restart[t]} {cfg.time_units[0]}")
+ if name_lower == "grid" and cfg.hide_map_elements[3] == 0 and cfg.title[k] == "0":
+ axis.set_title(
+ f"Grid = [{nx},{ny},{nz}], "
+ + f"Total no. active cells = {np.max(actind)+1}"
+ )
+ if cfg.title[k] != "0" and cfg.hide_map_elements[3] == 0 and not cfg.gif:
+ axis.set_title(cfg.title[k])
+ if cfg.slice[n_s][2][0] == -2 and not axis.yaxis_inverted():
+ axis.invert_yaxis()
+ if len(cfg.xlim[n]) > 1:
+ axis.set_xlim([float(cfg.xlim[n][0][1:]), float(cfg.xlim[n][1][:-1])])
+ xlabels = np.linspace(
+ float(cfg.xlim[n][0][1:]) * cfg.xscale,
+ float(cfg.xlim[n][1][:-1]) * cfg.xscale,
+ int(cfg.xtick_count[n]),
+ )
+ else:
+ xlabels = np.linspace(
+ np.min(xc) * cfg.xscale,
+ np.max(xc) * cfg.xscale,
+ int(cfg.xtick_count[n]),
+ )
+ _set_axis_ticks(
+ axis, "x", xlabels, cfg.xscale, cfg.xformat[n], cfg.hide_map_elements[1]
+ )
+ if len(cfg.ylim[n]) > 1:
+ axis.set_ylim([float(cfg.ylim[n][0][1:]), float(cfg.ylim[n][1][:-1])])
+ ylabels = np.linspace(
+ float(cfg.ylim[n][0][1:]) * cfg.yscale,
+ float(cfg.ylim[n][1][:-1]) * cfg.yscale,
+ int(cfg.ytick_count[n]),
+ )
+ else:
+ ylabels = np.linspace(
+ np.min(yc) * cfg.yscale,
+ np.max(yc) * cfg.yscale,
+ int(cfg.ytick_count[n]),
+ )
+ _set_axis_ticks(
+ axis, "y", ylabels, cfg.yscale, cfg.yformat[n], cfg.hide_map_elements[0]
+ )
+
+
+def _formatted_ticks(
+ values: NDArray,
+ scale: float,
+ value_format: str,
+) -> tuple[list[float], list[str]]:
+ """Format scaled tick locations and labels.
+
+ Parameters
+ ----------
+ values : np.ndarray
+ Unscaled tick values.
+ scale : float
+ Coordinate scale factor.
+ value_format : str
+ Python format specification.
+
+ Returns
+ -------
+ tuple[list[float], list[str]]
+ Scaled tick locations and formatted labels.
+
+ """
+ labels = [format(value, value_format) for value in values]
+ ticks = [float(label) / scale for label in labels]
+ return ticks, labels
+
+
+def _set_axis_ticks(
+ axis: Any,
+ axis_name: str,
+ labels: NDArray,
+ scale: float,
+ value_format: str,
+ remove_axis: int,
+) -> None:
+ """Set formatted ticks on one coordinate axis.
+
+ Parameters
+ ----------
+ axis : matplotlib.axes.Axes or np.ndarray
+ Map axes.
+ axis_name : {"x", "y"}
+ Coordinate axis to update.
+ labels : np.ndarray
+ Tick values before scaling.
+ scale : float
+ Coordinate scale factor.
+ value_format : str
+ Python format specification.
+ remove_axis : int
+ Nonzero when the selected axis is hidden.
+
+ """
+ if axis_name == "x":
+ if value_format and remove_axis == 0:
+ ticks, ticklabels = _formatted_ticks(labels, scale, value_format)
+ axis.set_xticks(ticks)
+ axis.set_xticklabels(ticklabels)
+ elif remove_axis == 0:
+ axis.set_xticks(labels / scale)
+ if scale != 1:
+ axis.set_xticklabels(labels)
+ else:
+ if value_format and remove_axis == 0:
+ ticks, ticklabels = _formatted_ticks(labels, scale, value_format)
+ axis.set_yticks(ticks)
+ axis.set_yticklabels(ticklabels)
+ elif remove_axis == 0:
+ axis.set_yticks(labels / scale)
+ if scale != 1:
+ axis.set_yticklabels(labels)
+
+
+def _draw_map(
+ deck: str,
+ fig: Figure,
+ axes: Any,
+ original_loc: list[Any],
+ cb: list[Any],
+ cmin: list[float],
+ cmax: list[float],
+ maska: list[Any],
+ diffa: list[NDArray],
+ named: str,
+ deckd: str,
+ slice_title: str,
+ slice_name: str,
+ cfg: PlopmConfig,
+ generated_files: list[str],
+ data: SimData,
+ t: int,
+ n: int,
+ k: int,
+ xc: NDArray,
+ yc: NDArray,
+ sub1: int,
+ mx: int,
+ my: int,
+ xname: str,
+ yname: str,
+) -> None:
+ """Draw and optionally save one spatial map.
+
+ Parameters
+ ----------
+ deck : str
+ Simulation-case stem.
+ fig : matplotlib.figure.Figure
+ Figure receiving the map.
+ axes : matplotlib.axes.Axes or np.ndarray
+ Target axes.
+ original_loc, cb : list
+ Original axes locators and active colorbars.
+ cmin, cmax : list[float]
+ Configured color limits.
+ maska, diffa : list
+ Mapped masks and cached difference arrays.
+ named, deckd : str
+ Display names for the primary and difference cases.
+ slice_title, slice_name : str
+ Human-readable slice descriptions.
+ cfg : PlopmConfig
+ Initialized map configuration.
+ generated_files : list[str]
+ Generated filenames updated when a PNG is saved.
+ data : SimData
+ Loaded simulation data.
+ t, n, k : int
+ Restart-step, variable, and subplot indices.
+ xc, yc : np.ndarray
+ Coordinate meshes.
+ sub1 : int
+ Number of subplot columns.
+ mx, my : int
+ Mapped grid dimensions.
+ xname, yname : str
+ Coordinate-axis names.
+
+ """
+ var = cfg.variables[n]
+ unit, values = read_quantity(
+ deck,
+ data,
+ var,
+ data.steps[t],
+ float(cfg.scale_factor[n]),
+ cfg.mass_vars,
+ cfg.mass_vars + cfg.mass_fracs,
+ cfg.caprock_vars,
+ cfg.stress_coefficient,
+ cfg.filters[k],
+ cfg.gif,
+ cfg.min_threshold[n],
+ cfg.max_threshold[n],
+ cfg.csv_columns[k],
+ )
+ n_s, feature_id, features = 0, 1, None
+ labels: list[str] = []
+ if cfg.subplot_grid[0] and len(cfg.cases[0]) > 1:
+ n_s = k
+ if cfg.csv_columns[k][0]:
+ quaa = values
+ else:
+ if cfg.variables[0] == "wells":
+ features, labels = get_wells(cfg, k)
+ elif cfg.variables[0] == "faults":
+ features, labels = get_faults(cfg, k)
+ feature_id = len(labels) + 1
+ if cfg.slice[n_s][0][0] != -2:
+ quaa = map_yz(cfg, data, var, values, k, mx, my, features, feature_id)
+ elif cfg.slice[n_s][1][0] != -2:
+ quaa = map_xz(cfg, data, var, values, k, mx, my, features, feature_id)
+ else:
+ quaa = map_xy(cfg, data, var, values, k, mx, my, features, feature_id)
+ if cfg.difference_input:
+ quaa -= diffa[t]
+ if cfg.mask_variable:
+ mask = maska[k]
+ maxv = np.nanmax(mask)
+ mask_condition = quaa < cfg.mask_threshold
+ quaa[mask_condition] = -cmax[n] * (maxv - mask[mask_condition]) / (maxv - 1)
+ if cfg.csv:
+ text = [f"{val}\n" for val in quaa if not np.isnan(val)]
+ name = _clean_name(f"{named}_{var}_{slice_name}_t{data.steps[t]}")
+ if cfg.filename[n]:
+ name = cfg.filename[n]
+ filename = f"{name}.csv"
+ with open(
+ f"{cfg.output_dir}/{filename}",
+ "w",
+ encoding="utf8",
+ ) as file:
+ file.write("".join(text))
+ generated_files.append(filename)
+ return
+ if var in cfg.mass_vars and cfg.difference_input:
+ extinf = np.nansum(np.abs(quaa))
+ elif var in cfg.mass_vars:
+ extinf = np.sum(quaa[~np.isnan(quaa)])
+ elif cfg.difference_input:
+ extinf = np.nansum(np.abs(quaa))
+ else:
+ extinf = np.empty(0)
+ ntick = 3
+ ncolor = var + " " + unit
+ defcol, temp, cmap = True, "tab20", matplotlib.colormaps.get_cmap("tab20")
+ if cfg.colormaps[n] in plt.colormaps():
+ defcol = False
+ cmap = matplotlib.colormaps.get_cmap(cfg.colormaps[n])
+ temp = cfg.colormaps[n]
+ if var not in ("wells", "grid", "faults"):
+ valid_maps = quaa[~np.isnan(quaa)]
+ if (
+ len(cfg.cases[0]) > 1
+ and cfg.subplot_grid[0]
+ or len(cfg.variables) > 1
+ and cfg.subplot_grid[0]
+ or len(data.steps) > 1
+ and cfg.subplot_grid[0]
+ and len(cfg.cases[0]) == 1
+ or cfg.gif
+ and not cfg.subplot_grid[0]
+ or int(cfg.color_log[n]) == 1
+ ):
+ minc = cmin[n]
+ maxc = cmax[n]
+ elif not cfg.global_range and valid_maps.size > 0:
+ minc = np.min(valid_maps)
+ maxc = np.max(valid_maps)
+ elif valid_maps.size > 0:
+ values = np.asarray(values)
+ valid_quan = values[~np.isnan(values)]
+ if valid_quan.size > 0:
+ minc = np.min(valid_quan)
+ maxc = np.max(valid_quan)
+ else:
+ minc = 0
+ maxc = 0
+ else:
+ minc = 0
+ maxc = 0
+ if cfg.clim[n][0]:
+ minc = float(cfg.clim[n][0][1:])
+ maxc = float(cfg.clim[n][1][:-1])
+ elif cfg.difference_input and int(cfg.color_log[n]) == 0:
+ minmax = max(abs(maxc), abs(minc))
+ minc = -minmax
+ maxc = minmax
+ if maxc == minc:
+ ntick = 1
+ elif (
+ "num" in var
+ and (cfg.colormaps[n] in cfg.disc_colormaps or defcol)
+ and cfg.discrete
+ and (minc.is_integer() and maxc.is_integer())
+ ):
+ ntick = int(maxc - minc + 1)
+ if cfg.mask_variable:
+ minc = -maxc
+ elif var in ["faults", "wells"]:
+ minc = 1
+ maxc = feature_id
+ else:
+ minc = 1
+ maxc = 1
+ nlc = ntick
+ if cfg.colorbar_tick_count[n] and ntick > 1:
+ ntick = int(cfg.colorbar_tick_count[n])
+ if cfg.colorbar_label:
+ ncolor = cfg.colorbar_label
+ shc = 0.0
+ if abs(minc) < sys.float_info.epsilon:
+ minc = 0
+ if ("num" in var and temp in cfg.disc_colormaps and cfg.discrete) or (
+ defcol and temp != "nipy_spectral"
+ ):
+ if maxc == minc:
+ shc = 2.0
+ from_list = matplotlib.colors.LinearSegmentedColormap.from_list
+ cmap = from_list(
+ "custom",
+ matplotlib.colormaps[temp](range(int(minc), int(minc) + nlc + int(shc))),
+ nlc,
+ )
+ if ntick == 2:
+ shc = (maxc - minc) / 2.0
+ elif minc == 0 and "num" not in var and var != "mpi_rank" or cfg.mask_variable:
+ shc = 0
+ else:
+ shc = 0.5
+ if defcol:
+ temp0 = []
+ for values in cfg.colormaps[n].split(" "):
+ if values[0] == "#":
+ temp0.append(values)
+ else:
+ temp0.append([])
+ for color in values.split(";"):
+ if color.isnumeric():
+ temp0[-1].append(float(color) / 255.0)
+ else:
+ plopm_error(
+ f"Color given in {cli_error_value(f'-c {cfg.colormaps[n]}')} not found."
+ )
+ cmap = colors.ListedColormap(temp0)
+ if cfg.inactive_color != "w":
+ cmap = cmap.with_extremes(bad=cfg.inactive_color)
+ axis = axes.flat[k]
+ if len(cfg.grid_edges) > 1:
+ if var == "grid":
+ imag = axis.pcolormesh(
+ xc,
+ yc,
+ quaa.reshape(my, mx),
+ facecolors="none",
+ edgecolors=cfg.grid_edges[0],
+ lw=float(cfg.grid_edges[1]),
+ )
+ elif int(cfg.color_log[n]) == 0:
+ imag = axis.pcolormesh(
+ xc,
+ yc,
+ quaa.reshape(my, mx),
+ shading="flat",
+ cmap=cmap,
+ edgecolors=cfg.grid_edges[0],
+ lw=float(cfg.grid_edges[1]),
+ )
+ else:
+ imag = axis.pcolormesh(
+ xc,
+ yc,
+ quaa.reshape(my, mx),
+ shading="flat",
+ cmap=cmap,
+ norm=colors.LogNorm(vmin=minc, vmax=maxc),
+ edgecolors=cfg.grid_edges[0],
+ lw=float(cfg.grid_edges[1]),
+ )
+ else:
+ if var == "grid":
+ imag = axis.pcolormesh(
+ xc,
+ yc,
+ quaa.reshape(my, mx),
+ facecolors="none",
+ edgecolors="black",
+ lw=0.001,
+ )
+ elif int(cfg.color_log[n]) == 0:
+ imag = axis.pcolormesh(
+ xc,
+ yc,
+ quaa.reshape(my, mx),
+ shading="flat",
+ cmap=cmap,
+ )
+ else:
+ imag = axis.pcolormesh(
+ xc,
+ yc,
+ quaa.reshape(my, mx),
+ shading="flat",
+ cmap=cmap,
+ norm=colors.LogNorm(vmin=minc, vmax=maxc),
+ )
+ if cfg.subplot_grid[0] and cfg.gif and len(cfg.variables) > 1 and cb[k] != "":
+ axes, cb = _remove_colorbar(axes, original_loc, cb, k)
+ if cfg.subplot_grid[0] and cfg.gif and len(cfg.cases[0]) > 1 and cb[k] != "":
+ axes, cb = _remove_colorbar(axes, original_loc, cb, k)
+ if (
+ not cfg.subplot_grid[0]
+ and cb[k] != ""
+ and cfg.gif
+ and cfg.hide_map_elements[2] == 0
+ ):
+ axes, cb = _remove_colorbar(axes, original_loc, cb, k)
+ divider = make_axes_locatable(axis)
+ if cfg.mask_variable:
+ vect = np.linspace(
+ 0,
+ maxc,
+ ntick,
+ endpoint=True,
+ )
+ else:
+ vect = np.linspace(
+ minc,
+ maxc,
+ ntick,
+ endpoint=True,
+ )
+ frmt = "{:" + cfg.cb_formats[n] + "}"
+
+ def formatter(value: float, _: Any) -> str:
+ """Format a colorbar value.
+
+ Parameters
+ ----------
+ value : float
+ Colorbar value.
+ _ : Any
+ Unused Matplotlib tick position.
+
+ Returns
+ -------
+ str
+ Formatted colorbar label.
+
+ """
+ return frmt.format(value)
+
+ if not cfg.mask_variable:
+ if int(cfg.color_log[n]) == 1:
+ pass
+ else:
+ for i, val in enumerate(vect):
+ if abs(float(frmt.format(val))) == 0:
+ vect[i] = 0
+ if i == 0:
+ minc = 0
+ if var not in ("wells", "grid", "faults"):
+ if int(cfg.color_log[n]) == 0:
+ if len(data.steps) > 1 and cfg.subplot_grid[0] and len(cfg.cases[0]) == 1:
+ if cfg.colorbar_position[0] != -1:
+ cb[0] = fig.colorbar(
+ imag,
+ cax=fig.add_axes(cfg.colorbar_position),
+ ticks=vect,
+ label=ncolor,
+ format=(
+ mticker.FixedFormatter(cfg.colorbar_ticks[n])
+ if cfg.colorbar_ticks[n]
+ else formatter
+ ),
+ shrink=0.2,
+ location="top",
+ )
+ elif not cfg.subplot_grid[0] or len(cfg.cases[0]) == 1:
+ cb[k] = fig.colorbar(
+ imag,
+ cax=divider.append_axes("right", size="2%", pad=0.05),
+ orientation="vertical",
+ ticks=vect,
+ label=ncolor,
+ format=(
+ mticker.FixedFormatter(cfg.colorbar_ticks[n])
+ if cfg.colorbar_ticks[n]
+ else formatter
+ ),
+ )
+ elif k == 0 and cfg.colorbar_position[0] != -1:
+ cb[0] = fig.colorbar(
+ imag,
+ cax=fig.add_axes(cfg.colorbar_position),
+ ticks=vect,
+ label=ncolor,
+ format=(
+ mticker.FixedFormatter(cfg.colorbar_ticks[n])
+ if cfg.colorbar_ticks[n]
+ else formatter
+ ),
+ shrink=0.2,
+ location="top",
+ )
+ else:
+ if cfg.color_log_ticks:
+
+ class LogTickFormatter(LogFormatter):
+ def set_locs(self, locs: Any | None = None) -> None:
+ """Set logarithmic colorbar sublabels from the configuration.
+
+ Parameters
+ ----------
+ locs : Any, optional
+ Tick locations supplied by Matplotlib.
+
+ """
+ self._sublabels = set(cfg.color_log_ticks)
+
+ if cfg.subplot_grid[0]:
+ if cfg.colorbar_position[0] != -1:
+ if cfg.color_log_ticks:
+ cb[k] = fig.colorbar(
+ imag,
+ cax=fig.add_axes(cfg.colorbar_position),
+ label=ncolor,
+ shrink=0.2,
+ location="top",
+ ticks=cfg.color_log_ticks,
+ format=LogTickFormatter(),
+ )
+ else:
+ cb[k] = fig.colorbar(
+ imag,
+ cax=fig.add_axes(cfg.colorbar_position),
+ label=ncolor,
+ shrink=0.2,
+ location="top",
+ )
+ else:
+ if cfg.color_log_ticks:
+ cb[k] = fig.colorbar(
+ imag,
+ cax=divider.append_axes("right", size="5%", pad=0.05),
+ orientation="vertical",
+ label=ncolor,
+ ticks=cfg.color_log_ticks,
+ format=LogTickFormatter(),
+ )
+ else:
+ cb[k] = fig.colorbar(
+ imag,
+ cax=divider.append_axes("right", size="5%", pad=0.05),
+ orientation="vertical",
+ label=ncolor,
+ )
+ else:
+ _add_map_overlay(fig, cfg, imag, divider, vect, n, var, features, labels)
+ imag.set_clim(
+ minc - shc,
+ maxc + shc,
+ )
+ _set_axis(
+ fig,
+ axes,
+ cfg,
+ data,
+ var,
+ n,
+ t,
+ k,
+ n_s,
+ unit,
+ xc,
+ yc,
+ extinf,
+ named,
+ deckd,
+ defcol,
+ slice_title,
+ feature_id,
+ )
+ if cfg.xlabel[n] and cfg.hide_map_elements[1] == 0:
+ axis.set_xlabel(cfg.xlabel[n])
+ elif (
+ cfg.hide_map_elements[1] == 0
+ and len(cfg.variables) == 1
+ and (k + sub1 >= len(cfg.cases[0]) or not cfg.subplot_grid[0])
+ ):
+ if len(data.steps) > 1 and cfg.subplot_grid[0] and len(cfg.cases[0]) == 1:
+ if k + sub1 >= len(data.steps):
+ axis.set_xlabel(f"{xname+cfg.xunit}")
+ else:
+ axis.set_xlabel(f"{xname+cfg.xunit}")
+ elif (
+ cfg.hide_map_elements[1] == 0
+ and len(cfg.cases[0]) == 1
+ and (k + sub1 >= len(cfg.variables) or not cfg.subplot_grid[0])
+ ) or (
+ cfg.hide_map_elements[1] == 0
+ and len(cfg.cases[0]) == len(cfg.variables)
+ and len(cfg.variables) > 1
+ and (k + sub1 >= len(cfg.variables) or not cfg.subplot_grid[0])
+ ):
+ axis.set_xlabel(f"{xname+cfg.xunit}")
+ if cfg.ylabel[n] and cfg.hide_map_elements[0] == 0:
+ axis.set_ylabel(cfg.ylabel[n])
+ elif cfg.hide_map_elements[0] == 0 and (k % sub1 == 0 or not cfg.subplot_grid[0]):
+ axis.set_ylabel(f"{yname+cfg.yunit}")
+ if cfg.hide_map_elements[2] == 1 and len(fig.axes) > 1:
+ fig.delaxes(fig.axes[1])
+ if (
+ cfg.hide_map_elements[1] == 1
+ or (
+ k + sub1 < len(cfg.cases[0])
+ and cfg.subplot_grid[0]
+ and len(cfg.variables) == 1
+ and cfg.remove_duplicate_labels
+ )
+ or cfg.hide_map_elements[1] == 1
+ or (
+ k + sub1 < len(cfg.variables)
+ and cfg.subplot_grid[0]
+ and len(cfg.cases[0]) == 1
+ and cfg.remove_duplicate_labels
+ )
+ or (
+ k + sub1 < len(data.steps)
+ and len(data.steps) > 1
+ and cfg.subplot_grid[0]
+ and len(cfg.cases[0]) == 1
+ and cfg.remove_duplicate_labels
+ )
+ ):
+ axis.tick_params(axis="x", which="both", bottom=False, labelbottom=False)
+ if cfg.hide_map_elements[0] == 1 or (
+ k % sub1 > 0 and cfg.subplot_grid[0] and cfg.remove_duplicate_labels == 1
+ ):
+ axis.tick_params(axis="y", which="both", left=False, labelleft=False)
+ axis.set_facecolor(cfg.fc)
+ if not cfg.gif:
+ if cfg.subplot_grid[0]:
+ if (
+ t == len(data.steps) - 1
+ and len(data.steps) > 1
+ or n == len(cfg.variables) - 1
+ and len(cfg.variables) > 1
+ ):
+ _save_map(
+ fig,
+ cfg,
+ data,
+ generated_files,
+ named,
+ var,
+ slice_name,
+ t,
+ n,
+ )
+ else:
+ if len(data.steps) == 1:
+ if k == max(len(cfg.variables) - 1, len(cfg.cases[0]) - 1):
+ _save_map(
+ fig,
+ cfg,
+ data,
+ generated_files,
+ named,
+ var,
+ slice_name,
+ t,
+ n,
+ )
+ elif (
+ len(cfg.cases[0]) == 1
+ or len(data.steps) > 1
+ and len(cfg.cases[0]) == len(data.steps)
+ ):
+ if t == len(data.steps) - 1:
+ _save_map(
+ fig,
+ cfg,
+ data,
+ generated_files,
+ named,
+ var,
+ slice_name,
+ t,
+ n,
+ )
+ else:
+ _save_map(
+ fig,
+ cfg,
+ data,
+ generated_files,
+ named,
+ var,
+ slice_name,
+ t,
+ n,
+ )
+ else:
+ save_index = t if cfg.rst_range else n
+ _save_map(
+ fig,
+ cfg,
+ data,
+ generated_files,
+ named,
+ var,
+ slice_name,
+ t,
+ save_index,
+ )
+ plt.close()
+
+
+def _clean_name(name: str) -> str:
+ """Convert a variable expression to a filename-safe stem.
+
+ Parameters
+ ----------
+ name : str
+ Variable expression or filename stem.
+
+ Returns
+ -------
+ str
+ Name with operators and spaces replaced.
+
+ """
+ name = name.replace(" / ", "_over_")
+ name = name.replace(" ", "")
+ return name
+
+
+def _save_map(
+ fig: Figure,
+ cfg: PlopmConfig,
+ data: SimData,
+ generated_files: list[str],
+ named: str,
+ var: str,
+ slice_name: str,
+ t: int,
+ save_index: int,
+) -> None:
+ """Save the current spatial map as a PNG file.
+
+ Parameters
+ ----------
+ fig : matplotlib.figure.Figure
+ Figure containing the map.
+ cfg : PlopmConfig
+ Output filename, directory, resolution, and face-color settings.
+ data : SimData
+ Simulation data containing the selected restart steps.
+ generated_files : list[str]
+ Generated filenames updated in place.
+ named : str
+ Case name used in the default filename.
+ var : str
+ Plotted variable name or expression.
+ slice_name : str
+ Slice description used in the default filename.
+ t : int
+ Index of the restart step being plotted.
+ save_index : int
+ Index used to select a custom filename.
+
+ """
+ fig.set_facecolor(cfg.fc)
+ name = _clean_name(f"{named}_{var}_{slice_name}_t{data.steps[t]}")
+ if save_index < len(cfg.filename) and cfg.filename[save_index]:
+ name = cfg.filename[save_index]
+
+ filename = f"{name}.png"
+ fig.savefig(
+ f"{cfg.output_dir}/{filename}",
+ bbox_inches="tight",
+ dpi=int(cfg.dpi[0]),
+ )
+ generated_files.append(filename)
+
+
+def _remove_colorbar(
+ axes: Any,
+ original_loc: list[Any],
+ cb: list[Any],
+ colorbar_index: int,
+) -> tuple[Any, list[Any]]:
+ """Remove a colorbar and restore its axes locator.
+
+ Parameters
+ ----------
+ axes : matplotlib.axes.Axes or np.ndarray
+ Map axes.
+ original_loc : list
+ Original axes locators.
+ cb : list
+ Active colorbar objects.
+ colorbar_index : int
+ Colorbar and axes index to restore.
+
+ Returns
+ -------
+ tuple
+ Updated axes and colorbar list.
+
+ """
+ if (
+ colorbar_index < len(cb)
+ and colorbar_index < len(original_loc)
+ and cb[colorbar_index] != ""
+ ):
+ cb[colorbar_index].remove()
+ axes.flat[colorbar_index].set_axes_locator(original_loc[colorbar_index])
+ cb[colorbar_index] = ""
+ return axes, cb
+
+
+def _add_map_overlay(
+ fig: Figure,
+ cfg: PlopmConfig,
+ imag: ScalarMappable,
+ divider: AxesDivider,
+ vect: NDArray,
+ n: int,
+ var: str,
+ features: list | None,
+ labels: list[str],
+) -> None:
+ """Add a categorical colorbar and labels for map features.
+
+ Parameters
+ ----------
+ fig : matplotlib.figure.Figure
+ Figure containing the map.
+ cfg : PlopmConfig
+ Initialized map configuration.
+ imag : matplotlib.cm.ScalarMappable
+ Mappable used to construct the colorbar.
+ divider : mpl_toolkits.axes_grid1.axes_divider.AxesDivider
+ Divider associated with the map axis.
+ vect : np.ndarray
+ Categorical colorbar tick values.
+ n : int
+ Variable index.
+ var : str
+ Categorical variable, such as ``"wells"`` or ``"faults"``.
+ features : list, optional
+ Feature locations grouped by label.
+ labels : list[str]
+ Feature names.
+
+ """
+ fig.colorbar(
+ imag,
+ cax=divider.append_axes("right", size="0%", pad=0.05),
+ orientation="vertical",
+ ticks=vect,
+ format=lambda x, _: "",
+ )
+ feature_id = len(labels) + 1
+ if var in ["faults", "wells"] and features is not None:
+ cmap = matplotlib.colormaps[cfg.colormaps[n]]
+ colour = cmap(np.linspace(0, 1, feature_id))
+ if feature_id < 70:
+ for label_index, label_name in enumerate(labels):
+ _add_label(features, label_index, label_name, colour)
+ else:
+ for label_index, label_name in zip(
+ [0, len(getattr(cfg, var)) - 1], [labels[0], labels[-1]]
+ ):
+ _add_label(features, label_index, label_name, colour)
+
+
+def _add_label(
+ features: list,
+ label_index: int,
+ label_name: str,
+ colour: NDArray,
+) -> None:
+ """Draw one feature label when the feature is present.
+
+ Parameters
+ ----------
+ features : list
+ Feature locations grouped by label.
+ label_index : int
+ Index of the feature group.
+ label_name : str
+ Text shown beside the categorical colorbar.
+ colour : np.ndarray
+ Colors assigned to feature groups.
+
+ """
+ item = features[label_index]
+ if any(item):
+ plt.text(
+ 0,
+ label_index + 1,
+ f"{label_name}",
+ c=colour[label_index],
+ fontweight="bold",
+ )
+
+
+
+
+
+# SPDX-FileCopyrightText: 2024-2026 NORCE Research AS
+# SPDX-License-Identifier: GPL-3.0
+# pylint: disable=W3301,R0912,R0913,R0914,R0915,R0917,E1102
+
+"""Create VTK files from OPM Flow simulation results.
+
+The module runs a minimal OPM Flow job when grid geometry is unavailable,
+populates VTU cell-data arrays for selected restart steps, and writes the PVD
+collection used to open the resulting time series.
+"""
+
+import os
+import shlex
+import shutil
+import sys
+from contextlib import nullcontext
+from subprocess import run
+
+import numpy as np
+from alive_progress import alive_bar
+from numpy.typing import NDArray
+
+from plopm.config.config import SimData
+from plopm.utils.readers import read_case, read_quantity
+from plopm.utils.terminal import cli_error_value, plopm_error, plopm_warning
+
+VTK_DTYPES = {
+ "Float64": np.float64,
+ "Float32": np.float32,
+ "Float16": np.float16,
+ "Int64": np.int64,
+ "UInt64": np.uint64,
+ "Int32": np.int32,
+ "UInt32": np.uint32,
+ "Int16": np.int16,
+ "UInt16": np.uint16,
+ "Int8": np.int8,
+ "UInt8": np.uint8,
+}
+
+
+
+[docs]
+def make_vtks(
+ flow: str,
+ names: list,
+ output: str,
+ save: list,
+ restart: list,
+ variables: list,
+ vtkformat_list: list,
+ vtknames: list,
+ gif: bool,
+ vtk: bool,
+ filters: list,
+ scales: list[str],
+ mass: list[str],
+ mass_all: list[str],
+ caprock: list[str],
+ stress: float,
+ filterss: list[str],
+) -> list:
+ """Create VTK time-series output for the configured cases.
+
+ A minimal OPM Flow run creates the grid-only VTU file when needed. Selected
+ properties are then read from INIT or UNRST output and written to one VTU
+ file per restart step.
+
+ Parameters
+ ----------
+ flow : str
+ Command used to run OPM Flow.
+ names : list
+ Simulation-case stems grouped by the CLI input.
+ output : str
+ Directory in which VTK files are written.
+ save : list
+ Optional output stems for each case.
+ restart : list
+ Restart report steps to export.
+ variables : list
+ Variables or expressions written as cell data.
+ vtkformat_list : list
+ VTK data type selected for each variable.
+ vtknames : list
+ Optional VTK array names for each variable.
+ gif, vtk : bool
+ Output-mode flags passed to the simulation readers.
+ filters : list
+ Property filters used while loading each case.
+ scales : list[str]
+ Scale factor applied to each variable.
+ mass, mass_all : list[str]
+ Mass variables and all supported mass-related variables.
+ caprock : list[str]
+ Supported caprock-integrity variables.
+ stress : float
+ Vertical stress coefficient used for caprock quantities.
+ filterss : list[str]
+ Filter expressions applied while reading exported quantities.
+
+ Returns
+ -------
+ list[str]
+ Names of the generated PVD collection files.
+
+ """
+ generated_files: list[str] = []
+
+ for k, case in enumerate(names[0]):
+ deck = case
+ dname = case.split("/")[-1]
+ grid_name = f"{dname}-GRID.vtu"
+ grid_path = os.path.join(output, grid_name)
+
+ if not os.path.isfile(f"{deck}.DATA"):
+ plopm_error(f"unable to find {cli_error_value(f'{deck}.DATA')}.")
+
+ if not os.path.isfile(grid_path):
+ cwd = os.getcwd()
+ output_abs = os.path.abspath(output)
+ dryrun_deck = ""
+ dryrun_folder = ""
+ dryrun_parent = cwd
+ try:
+ if len(case.split("/")) > 1:
+ os.chdir("/".join(case.split("/")[:-1]))
+ dryrun_parent = os.getcwd()
+ flags, thermal = _vtk_flags()
+ flow_command = shlex.split(flow)
+ dryrun_deck = f"{dname}_DRYRUN_{os.getpid()}.DATA"
+ dryrun_folder = f"plopm_{os.getpid()}"
+ shutil.copyfile(f"{dname}.DATA", dryrun_deck)
+ flags += " --enable-dry-run=1"
+ os.makedirs(dryrun_folder, exist_ok=True)
+ deck_rel = f"../{dryrun_deck}"
+ os.chdir(dryrun_folder)
+ if "SPE11B" in dname or "SPE11C" in dname:
+ run(
+ flow_command
+ + [deck_rel]
+ + shlex.split(flags)
+ + shlex.split(thermal),
+ check=False,
+ )
+ else:
+ run(flow_command + [deck_rel] + shlex.split(flags), check=False)
+ shutil.move(
+ f"{dname}_DRYRUN_{os.getpid()}-00000.vtu",
+ os.path.join(output_abs, grid_name),
+ )
+ finally:
+ os.chdir(dryrun_parent)
+ if dryrun_folder:
+ shutil.rmtree(dryrun_folder, ignore_errors=True)
+ if dryrun_deck and os.path.isfile(dryrun_deck):
+ os.remove(dryrun_deck)
+ os.chdir(cwd)
+
+ generated_files.append(grid_name)
+
+ data = read_case(case, gif, vtk, variables, restart, filters)
+ _write_vtk_data(
+ case,
+ data,
+ output,
+ dname,
+ save,
+ variables,
+ vtkformat_list,
+ vtknames,
+ k,
+ scales,
+ mass,
+ mass_all,
+ caprock,
+ stress,
+ filterss[k],
+ )
+
+ where = save[k] if save[k] else dname
+ generated_files.extend(
+ f"{where}-{int(restart_index):04d}.vtu" for restart_index in data.steps
+ )
+
+ _write_pvd(
+ save,
+ dname,
+ data.steps,
+ data.times,
+ output,
+ k,
+ )
+ generated_files.append(f"{where}.pvd")
+
+ return list(dict.fromkeys(generated_files))
+
+
+
+def _write_pvd(
+ save: list, dname: str, restart: list, tnrst: list, output: str, k: int
+) -> None:
+ """Write a PVD collection for a VTU time series.
+
+ Parameters
+ ----------
+ save : list
+ Optional output stems for each case.
+ dname : str
+ Default case name.
+ restart : list
+ Restart report steps included in the collection.
+ tnrst : list
+ Simulation times indexed by restart report step.
+ output : str
+ Output directory.
+ k : int
+ Case index used to select the output stem.
+
+ """
+ where = save[k] if save[k] else dname
+ pvd_lines = []
+ pvd_lines.append(
+ "<?xml version='1.0'?>\n"
+ + "<VTKFile type='Collection'\n"
+ + " version='0.1'\n"
+ + " byte_order='LittleEndian'\n"
+ + " compressor='vtkZLibDataCompressor'>\n"
+ + " <Collection>\n"
+ )
+ for i in restart:
+ pvd_lines.append(
+ f" <DataSet timestep='{tnrst[i]}' file='{where}-{int(i):04d}.vtu'/>\n"
+ )
+ pvd_lines.append(" </Collection>\n</VTKFile>")
+ with open(
+ f"{output}/{where}.pvd",
+ "w",
+ encoding="utf8",
+ ) as file:
+ file.write("".join(pvd_lines))
+
+
+def _warn_once(warning_keys: set, warning_key, message: str) -> None:
+ """Emit a warning once for a unique key.
+
+ Parameters
+ ----------
+ warning_keys : set
+ Keys for warnings already emitted.
+ warning_key
+ Hashable key identifying the warning condition.
+ message : str
+ Warning message.
+
+ """
+ if warning_key not in warning_keys:
+ plopm_warning(message)
+ warning_keys.add(warning_key)
+
+
+def _check_integer_conversion(
+ values: NDArray,
+ var: str,
+ vtkformat: str,
+ target_dtype: type,
+ warning_keys: set[tuple[str, str, str]],
+) -> None:
+ """Warn about unsafe conversion to an integer VTK type.
+
+ Warnings cover non-numeric or non-finite values, negative values converted
+ to unsigned integers, decimal truncation, and values outside the target
+ integer range.
+
+ Parameters
+ ----------
+ values : np.ndarray
+ Quantity values to inspect.
+ var : str
+ Variable name used in warning messages.
+ vtkformat : str
+ Requested VTK data type.
+ target_dtype : type
+ NumPy dtype used for conversion.
+ warning_keys : set[tuple[str, str, str]]
+ Keys for warnings already emitted.
+
+ """
+ try:
+ numeric_values = np.asarray(values, dtype=np.float64)
+ except TypeError:
+ _warn_once(
+ warning_keys,
+ (var.upper(), vtkformat, "non_numeric"),
+ f"{var.upper()} contains non-numeric values but is written as {vtkformat}.",
+ )
+ return
+ if not numeric_values.size:
+ return
+ finite_mask = np.isfinite(numeric_values)
+ finite_values = numeric_values[finite_mask]
+ if finite_values.size != numeric_values.size:
+ _warn_once(
+ warning_keys,
+ (var.upper(), vtkformat, "non_finite"),
+ f"{var.upper()} contains non-finite values but is written as {vtkformat}.",
+ )
+ if not finite_values.size:
+ return
+ dtype_info = np.iinfo(target_dtype)
+ min_val = finite_values.min()
+ max_val = finite_values.max()
+ if np.issubdtype(target_dtype, np.unsignedinteger) and min_val < 0:
+ _warn_once(
+ warning_keys,
+ (var.upper(), vtkformat, "negative_unsigned"),
+ f"{var.upper()} contains negative values but is written as {vtkformat}; "
+ "NumPy may wrap them.",
+ )
+ if np.any(finite_values != np.trunc(finite_values)):
+ _warn_once(
+ warning_keys,
+ (var.upper(), vtkformat, "float_truncation"),
+ f"{var.upper()} contains float values but is written as {vtkformat}; "
+ "NumPy will truncate decimals.",
+ )
+ if min_val < dtype_info.min or max_val > dtype_info.max:
+ _warn_once(
+ warning_keys,
+ (var.upper(), vtkformat, "out_of_range"),
+ f"{var.upper()} contains values outside {vtkformat} range [{dtype_info.min}, "
+ f"{dtype_info.max}]; NumPy may wrap or fail depending on version.",
+ )
+
+
+def _format_vtk_array(values: NDArray, target_dtype: type) -> str:
+ """Format values for an ASCII VTK DataArray.
+
+ Parameters
+ ----------
+ values : np.ndarray
+ Values to flatten and convert.
+ target_dtype : type
+ NumPy dtype used for the output values.
+
+ Returns
+ -------
+ str
+ Tab-indented values ready for insertion into a VTU file.
+
+ """
+ values = np.ravel(np.asarray(values, dtype=target_dtype))
+ if np.issubdtype(np.dtype(target_dtype), np.floating):
+ values = np.char.mod("%.8f", values)
+ values = np.char.rstrip(np.char.rstrip(values, "0"), ".")
+ values = np.where(values == "-0", "0", values)
+ else:
+ values = values.astype(str)
+ return "\t\t\t\t\t " + " ".join(values) + "\n\t\t\t\t\t</DataArray>"
+
+
+def _write_vtk_data(
+ case: str,
+ data: SimData,
+ output: str,
+ dname: str,
+ save: list,
+ variables: list,
+ vtkformat_list: list,
+ vtknames: list,
+ k: int,
+ scales: list[str],
+ mass: list[str],
+ mass_all: list[str],
+ caprock: list[str],
+ stress: float,
+ filterss: str,
+) -> None:
+ """Populate grid VTU files with simulation cell data.
+
+ Parameters
+ ----------
+ case : str
+ Simulation-case stem.
+ data : SimData
+ Loaded OPM simulation data.
+ output : str
+ Output directory.
+ dname : str
+ Default case name.
+ save : list
+ Optional output stems for each case.
+ variables : list
+ Variables or expressions written as cell data.
+ vtkformat_list : list
+ VTK data type selected for each variable.
+ vtknames : list
+ Optional VTK array names.
+ k : int
+ Case index used to select output settings.
+ scales : list[str]
+ Scale factor applied to each variable.
+ mass, mass_all : list[str]
+ Mass variables and all supported mass-related variables.
+ caprock : list[str]
+ Supported caprock-integrity variables.
+ stress : float
+ Vertical stress coefficient used for caprock quantities.
+ filterss : str
+ Filter expression applied while reading quantities.
+
+ """
+ restart = data.steps
+ vtk_lines = []
+ skip = False
+ warning_keys: set[tuple[str, str, str]] = set()
+ with open(f"{output}/{dname}-GRID.vtu", encoding="utf8") as file:
+ for line in file:
+ if skip and "CellData" in line:
+ skip = False
+ continue
+ if "CellData" in line:
+ skip = True
+ if not skip:
+ vtk_lines.append(line)
+ where = save[k] if save[k] else dname
+ show_progress = sys.stdout.isatty()
+ if show_progress:
+ bar_ctx = alive_bar(len(restart) * len(variables), bar="fish")
+ else:
+ bar_ctx = nullcontext()
+ with bar_ctx as bar_animation:
+ for i in restart:
+ cell_data = [
+ "\t\t\t\t<CellData Scalars='File created by https://github.com/cssr-tools/plopm'>",
+ ]
+ for n, var in enumerate(variables):
+ if show_progress:
+ bar_animation()
+ unit, values = read_quantity(
+ case,
+ data,
+ var,
+ i,
+ float(scales[n]),
+ mass,
+ mass_all,
+ caprock,
+ stress,
+ filterss,
+ False,
+ "",
+ "",
+ [False],
+ )
+
+ vtkformat = vtkformat_list[n]
+ target_dtype = VTK_DTYPES[vtkformat]
+ if np.issubdtype(target_dtype, np.integer):
+ _check_integer_conversion(
+ values, var, vtkformat, target_dtype, warning_keys
+ )
+ # VTK XML interoperability for Float16 is limited in many readers,
+ # so we emit Float32 in the DataArray type while preserving values.
+ if vtkformat == "Float16":
+ vtkformat = "Float32"
+ cell_data.append(
+ f"\n\t\t\t\t\t<DataArray type='{vtkformat}' Name="
+ + f"'{vtknames[n] if vtknames[n] else var+unit}' "
+ + "NumberOfComponents='1' format='ascii'>\n"
+ )
+ cell_data.append(_format_vtk_array(values, target_dtype))
+ cell_data.append("\n\t\t\t\t</CellData>\n")
+ with open(
+ f"{output}/{where}-{int(i):04d}.vtu",
+ "w",
+ encoding="utf8",
+ ) as file:
+ file.write("".join(vtk_lines[:4] + cell_data + vtk_lines[4:]))
+
+
+def _vtk_flags() -> tuple[str, str]:
+ """Build OPM Flow options for a minimal VTK run.
+
+ Returns
+ -------
+ tuple[str, str]
+ General VTK options and optional thermal-model options.
+
+ """
+ flags = (
+ " --enable-vtk-output=1 --enable-ecl-output=0 --output-mode=none"
+ + " --vtk-write-temperature=0 --vtk-write-densities=0 --vtk-write-mole-fractions=0 "
+ + "--vtk-write-relative-permeabilities=0 --vtk-write-pressures=0 "
+ + "--vtk-write-saturations=0 --vtk-write-porosity=0"
+ )
+ thermal = ""
+ return flags, thermal
+