diff --git a/CHANGELOG.md b/CHANGELOG.md index 7ce10b298..3e1db6735 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ This project adheres to [Semantic Versioning](https://semver.org/). ## [Unreleased] ### Added +- [#160](https://github.com/equinor/flownet/pull/160) Adds the possibility to extract regions from an existing model when the data source is a simulation model. For equil, the scheme key can be set to 'regions_from_sim' - [#157](https://github.com/equinor/flownet/pull/157) Adds a new 'time-weighted average open perforation location' perforation strategy called `time_avg_open_location`. - [#150](https://github.com/equinor/flownet/pull/150) Adds this changelog. - [#146](https://github.com/equinor/flownet/pull/146) Added about page to documentation with logo of industry and research institute partners. diff --git a/src/flownet/ahm/_run_ahm.py b/src/flownet/ahm/_run_ahm.py index 5beed4a1b..c0a4e4bb5 100644 --- a/src/flownet/ahm/_run_ahm.py +++ b/src/flownet/ahm/_run_ahm.py @@ -3,6 +3,7 @@ import numpy as np import pandas as pd +from scipy.stats import mode from scipy.optimize import minimize from configsuite import ConfigSuite @@ -22,6 +23,62 @@ from ..data import FlowData +def _from_regions_to_flow_tubes( + network: NetworkModel, + field_data: FlowData, + ti2ci: pd.DataFrame, + region_name: str, +) -> List[int]: + """ + The function loops through each cell in all flow tubes, and checks what region the + corresponding position (cell midpoint) in the data source simulation model has. If different + cells in one flow tube are located in different regions of the original model, the mode is used. + If flow tubes are entirely outside of the data source simulation grid, + the region closest to the midpoint of the flow tube is used. + + Args: + network: FlowNet network instance + field_data: FlowData class with information from simulation model data source + ti2ci: A dataframe with index equal to tube model index, and one column which equals cell indices. + name: The name of the region parameter + + Returns: + A list with values for 'name' region for each tube in the FlowNet model + """ + df_regions = [] + + xyz_mid = network.cell_midpoints + + for i in network.grid.model.unique(): + tube_regions = [] + for j in ti2ci[ti2ci.index == i].values: + ijk = field_data.grid.find_cell(xyz_mid[0][j], xyz_mid[1][j], xyz_mid[2][j]) + if ijk is not None and field_data.grid.active(ijk=ijk): + tube_regions.append(field_data.init(region_name)[ijk]) + if tube_regions != []: + df_regions.append(mode(tube_regions).mode.tolist()[0]) + else: + tube_midpoint = network.get_connection_midpoints(int(i)) + dist_to_cell = [] + for k in range(1, field_data.grid.get_num_active()): + cell_midpoint = field_data.grid.get_xyz(active_index=k) + dist_to_cell.append( + np.sqrt( + np.square(tube_midpoint[0] - cell_midpoint[0]) + + np.square(tube_midpoint[1] - cell_midpoint[1]) + + np.square(tube_midpoint[2] - cell_midpoint[2]) + ) + ) + df_regions.append( + field_data.init(region_name)[ + field_data.grid.get_ijk( + active_index=dist_to_cell.index(min(dist_to_cell)) + ) + ] + ) + return df_regions + + def _find_training_set_fraction( schedule: Schedule, config: ConfigSuite.snapshot ) -> float: @@ -202,12 +259,6 @@ def run_flownet_history_matching( df_satnum = pd.DataFrame( [1] * len(network.grid.model.unique()), columns=["SATNUM"] ) - else: - raise ValueError( - f"The relative permeability scheme " - f"'{config.model_parameters.relative_permeability.scheme}' is not valid.\n" - f"Valid options are 'global' or 'individual'." - ) # Create a pandas dataframe with all parameter definition for each individual tube relperm_dist_values = pd.DataFrame( @@ -224,7 +275,7 @@ def run_flownet_history_matching( key: relperm_dict[key] for key in relperm_dict if key != "scheme" } - for i in df_satnum["SATNUM"].unique(): + for i in np.sort(df_satnum["SATNUM"].unique()): info = [ relperm_parameters.keys(), [relperm_parameters[key].min for key in relperm_parameters], @@ -250,37 +301,62 @@ def run_flownet_history_matching( df_eqlnum = pd.DataFrame( range(1, len(network.grid.model.unique()) + 1), columns=["EQLNUM"] ) + elif config.model_parameters.equil.scheme == "regions_from_sim": + df_eqlnum = pd.DataFrame( + _from_regions_to_flow_tubes(network, field_data, ti2ci, "EQLNUM"), + columns=["EQLNUM"], + ) elif config.model_parameters.equil.scheme == "global": df_eqlnum = pd.DataFrame( [1] * len(network.grid.model.unique()), columns=["EQLNUM"] ) - else: - raise ValueError( - f"The equilibration scheme " - f"'{config.model_parameters.relative_permeability.scheme}' is not valid.\n" - f"Valid options are 'global' or 'individual'." - ) # Create a pandas dataframe with all parameter definition for each individual tube equil_dist_values = pd.DataFrame( columns=["parameter", "minimum", "maximum", "loguniform", "eqlnum"] ) - equil_config = config.model_parameters.equil - for i in df_eqlnum["EQLNUM"].unique(): + defined_eqlnum_regions = [] + datum_depths = [] + if config.model_parameters.equil.scheme == "regions_from_sim": + equil_config_eqlnum = config.model_parameters.equil.regions + for reg in equil_config_eqlnum: + defined_eqlnum_regions.append(reg.id) + else: + equil_config_eqlnum = [config.model_parameters.equil.regions[0]] + defined_eqlnum_regions.append(None) + + for i in np.sort(df_eqlnum["EQLNUM"].unique()): + if i in defined_eqlnum_regions: + idx = defined_eqlnum_regions.index(i) + else: + idx = defined_eqlnum_regions.index(None) + datum_depths.append(equil_config_eqlnum[idx].datum_depth) info = [ ["datum_pressure", "owc_depth", "gwc_depth", "goc_depth"], [ - equil_config.datum_pressure.min, - None if equil_config.owc_depth is None else equil_config.owc_depth.min, - None if equil_config.gwc_depth is None else equil_config.gwc_depth.min, - None if equil_config.goc_depth is None else equil_config.goc_depth.min, + equil_config_eqlnum[idx].datum_pressure.min, + None + if equil_config_eqlnum[idx].owc_depth is None + else equil_config_eqlnum[idx].owc_depth.min, + None + if equil_config_eqlnum[idx].gwc_depth is None + else equil_config_eqlnum[idx].gwc_depth.min, + None + if equil_config_eqlnum[idx].goc_depth is None + else equil_config_eqlnum[idx].goc_depth.min, ], [ - equil_config.datum_pressure.max, - None if equil_config.owc_depth is None else equil_config.owc_depth.max, - None if equil_config.gwc_depth is None else equil_config.gwc_depth.max, - None if equil_config.goc_depth is None else equil_config.goc_depth.max, + equil_config_eqlnum[idx].datum_pressure.max, + None + if equil_config_eqlnum[idx].owc_depth is None + else equil_config_eqlnum[idx].owc_depth.max, + None + if equil_config_eqlnum[idx].gwc_depth is None + else equil_config_eqlnum[idx].gwc_depth.max, + None + if equil_config_eqlnum[idx].goc_depth is None + else equil_config_eqlnum[idx].goc_depth.max, ], [False] * 4, [i] * 4, @@ -361,6 +437,8 @@ def run_flownet_history_matching( # ****************************************************************************** + datum_depths = list(datum_depths) + parameters = [ PorvPoroTrans(porv_poro_trans_dist_values, ti2ci, network), RelativePermeability( @@ -375,7 +453,7 @@ def run_flownet_history_matching( network, ti2ci, df_eqlnum, - equil_config.datum_depth, + datum_depths, config.flownet.pvt.rsvd, ), ] diff --git a/src/flownet/config_parser/_config_parser.py b/src/flownet/config_parser/_config_parser.py index f22b2ca76..7ca278608 100644 --- a/src/flownet/config_parser/_config_parser.py +++ b/src/flownet/config_parser/_config_parser.py @@ -6,6 +6,8 @@ import configsuite from configsuite import types, MetaKeys as MK, ConfigSuite +from ..data.from_flow import FlowData + # Small workaround while waiting for https://github.com/equinor/configsuite/pull/157 # to be merged and released in upstream ConfigSuite: @@ -26,6 +28,23 @@ def create_schema(config_folder: Optional[pathlib.Path] = None) -> Dict: """ + @configsuite.transformation_msg("Convert 'None' to None") + def _none_to_none( + input_data: Union[str, int, float, None] + ) -> Union[str, int, float, None]: + """ + Converts "None" to None + Args: + input_data: + + Returns: + The input_data. If the input is "None" it is converted to None (str to None) + """ + if input_data == "None": + return None + + return input_data + @configsuite.transformation_msg("Convert string to lower case") def _to_lower(input_data: Union[List[str], str]) -> Union[List[str], str]: if isinstance(input_data, str): @@ -652,58 +671,74 @@ def _to_abs_path(path: Optional[str]) -> str: "scheme": { MK.Type: types.String, MK.Description: "Number of degrees of freedom. Either " - "'global' (one for the whole model) or 'individual' " - "(one per tube)", + "'global' (the whole model treated as one EQLNUM region)," + "'individual' (each tube treated as an EQLNUM region) or" + "'regions_from_sim' (EQLNUM regions extracted from input simulation", MK.Default: "global", MK.Transformation: _to_lower, }, - "datum_depth": { - MK.Type: types.Number, - MK.AllowNone: True, - }, - "datum_pressure": { - MK.Type: types.NamedDict, + "regions": { + MK.Type: types.List, MK.Content: { - "min": {MK.Type: types.Number}, - "max": {MK.Type: types.Number}, - }, - }, - "owc_depth": { - MK.Type: types.NamedDict, - MK.Content: { - "min": { - MK.Type: types.Number, - MK.AllowNone: True, - }, - "max": { - MK.Type: types.Number, - MK.AllowNone: True, - }, - }, - }, - "gwc_depth": { - MK.Type: types.NamedDict, - MK.Content: { - "min": { - MK.Type: types.Number, - MK.AllowNone: True, - }, - "max": { - MK.Type: types.Number, - MK.AllowNone: True, - }, - }, - }, - "goc_depth": { - MK.Type: types.NamedDict, - MK.Content: { - "min": { - MK.Type: types.Number, - MK.AllowNone: True, - }, - "max": { - MK.Type: types.Number, - MK.AllowNone: True, + MK.Item: { + MK.Type: types.NamedDict, + MK.Content: { + "id": { + MK.Type: types.Number, + MK.AllowNone: True, + MK.Transformation: _none_to_none, + }, + "datum_depth": { + MK.Type: types.Number, + MK.AllowNone: True, + }, + "datum_pressure": { + MK.Type: types.NamedDict, + MK.Content: { + "min": {MK.Type: types.Number}, + "max": {MK.Type: types.Number}, + }, + }, + "owc_depth": { + MK.Type: types.NamedDict, + MK.Content: { + "min": { + MK.Type: types.Number, + MK.AllowNone: True, + }, + "max": { + MK.Type: types.Number, + MK.AllowNone: True, + }, + }, + }, + "gwc_depth": { + MK.Type: types.NamedDict, + MK.Content: { + "min": { + MK.Type: types.Number, + MK.AllowNone: True, + }, + "max": { + MK.Type: types.Number, + MK.AllowNone: True, + }, + }, + }, + "goc_depth": { + MK.Type: types.NamedDict, + MK.Content: { + "min": { + MK.Type: types.Number, + MK.AllowNone: True, + }, + "max": { + MK.Type: types.Number, + MK.AllowNone: True, + }, + }, + }, + }, }, }, }, @@ -789,7 +824,7 @@ def parse_config(configuration_file: pathlib.Path) -> ConfigSuite.snapshot: Parsed config, where values can be extracted like e.g. 'config.ert.queue.system'. """ - # pylint: disable=too-many-branches + # pylint: disable=too-many-branches, too-many-statements, too-many-lines input_config = yaml.safe_load(configuration_file.read_text()) @@ -808,6 +843,54 @@ def parse_config(configuration_file: pathlib.Path) -> ConfigSuite.snapshot: config = suite.snapshot req_relp_parameters: List[str] = [] + if ( + config.model_parameters.equil.scheme != "regions_from_sim" + and config.model_parameters.equil.scheme != "individual" + and config.model_parameters.equil.scheme != "global" + ): + raise ValueError( + f"The relative permeability scheme " + f"'{config.model_parameters.equil.scheme}' is not valid.\n" + f"Valid options are 'global', 'regions_from_sim' or 'individual'." + ) + + if config.model_parameters.equil.scheme == "regions_from_sim": + if config.flownet.data_source.simulation.input_case is None: + raise ValueError( + "Input simulation case is not defined - EQLNUM regions can not be extracted" + ) + field_data = FlowData(config.flownet.data_source.simulation.input_case) + unique_regions = field_data.get_unique_regions("EQLNUM") + default_exists = False + defined_regions = [] + for reg in config.model_parameters.equil.regions: + if reg.id is None: + default_exists = True + else: + if reg.id in defined_regions: + raise ValueError(f"EQLNUM region {reg.id} defined multiple times") + defined_regions.append(reg.id) + + if reg.id not in unique_regions and reg.id is not None: + raise ValueError( + f"EQLNUM regions {reg.id} is not found in the input simulation case" + ) + + if set(defined_regions) != set(unique_regions): + print( + "Values not defined for all EQLNUM regions. Default values will be used if defined." + ) + if not default_exists: + raise ValueError("Default values for EQLNUM regions not defined") + + if ( + config.model_parameters.equil.scheme != "regions_from_sim" + and config.model_parameters.equil.regions[0].id is not None + ): + raise ValueError( + "Id for first equilibrium region parameter should not be set, or set to 'None'\n" + "when using the 'global' or 'individual' options" + ) for phase in config.flownet.phases: if phase not in ["oil", "gas", "water", "disgas", "vapoil"]: @@ -838,15 +921,16 @@ def parse_config(configuration_file: pathlib.Path) -> ConfigSuite.snapshot: "krwend", "krowend", ] - if ( - config.model_parameters.equil.owc_depth.min is None - or config.model_parameters.equil.owc_depth.max is None - or config.model_parameters.equil.owc_depth.max - < config.model_parameters.equil.owc_depth.min - ): - raise ValueError( - "Ambiguous configuration input: OWC not properly specified. Min or max missing, or max < min." - ) + for reg in config.model_parameters.equil.regions: + if ( + reg.owc_depth.min is None + or reg.owc_depth.max is None + or reg.owc_depth.max < reg.owc_depth.min + ): + raise ValueError( + "Ambiguous configuration input: OWC not properly specified. Min or max missing, or max < min." + ) + if {"oil", "gas"}.issubset(config.flownet.phases): req_relp_parameters = req_relp_parameters + [ "scheme", @@ -859,15 +943,15 @@ def parse_config(configuration_file: pathlib.Path) -> ConfigSuite.snapshot: "krgend", "krogend", ] - if ( - config.model_parameters.equil.goc_depth.min is None - or config.model_parameters.equil.goc_depth.max is None - or config.model_parameters.equil.goc_depth.max - < config.model_parameters.equil.goc_depth.min - ): - raise ValueError( - "Ambiguous configuration input: GOC not properly specified. Min or max missing, or max < min." - ) + for reg in config.model_parameters.equil.regions: + if ( + reg.goc_depth.min is None + or reg.goc_depth.max is None + or reg.goc_depth.max < reg.goc_depth.min + ): + raise ValueError( + "Ambiguous configuration input: GOC not properly specified. Min or max missing, or max < min." + ) for parameter in set(req_relp_parameters): if parameter == "scheme": diff --git a/src/flownet/data/from_flow.py b/src/flownet/data/from_flow.py index 38312d587..a71ee6887 100644 --- a/src/flownet/data/from_flow.py +++ b/src/flownet/data/from_flow.py @@ -7,7 +7,7 @@ import pandas as pd from ecl.well import WellInfo from ecl.grid import EclGrid -from ecl.eclfile import EclFile +from ecl.eclfile import EclFile, EclInitFile from ecl.summary import EclSum from ecl2df import faults from ecl2df.eclfiles import EclFiles @@ -37,6 +37,7 @@ def __init__( self._eclsum = EclSum(str(self._input_case)) self._grid = EclGrid(str(self._input_case.with_suffix(".EGRID"))) self._restart = EclFile(str(self._input_case.with_suffix(".UNRST"))) + self._init = EclInitFile(self._grid, str(self._input_case.with_suffix(".INIT"))) self._wells = WellInfo( self._grid, rst_file=self._restart, load_segment_information=True ) @@ -324,6 +325,14 @@ def _grid_cell_bounding_boxes(self) -> np.ndarray: def _get_start_date(self): return self._eclsum.start_date + def init(self, name: str) -> np.ndarray: + """array with 'name' regions""" + return self._init[name][0] + + def get_unique_regions(self, name: str) -> np.ndarray: + """array with unique 'name' regions""" + return np.unique(self._init[name][0]) + @property def grid_cell_bounding_boxes(self) -> np.ndarray: """Boundingboxes for all gridcells""" @@ -343,3 +352,8 @@ def production(self) -> pd.DataFrame: def coordinates(self) -> pd.DataFrame: """dataframe with all coordinates""" return self._coordinates() + + @property + def grid(self) -> EclGrid: + """the simulation grid with properties""" + return self._grid diff --git a/src/flownet/network_model/_network_model.py b/src/flownet/network_model/_network_model.py index 8a445d043..1f1e61d0d 100644 --- a/src/flownet/network_model/_network_model.py +++ b/src/flownet/network_model/_network_model.py @@ -1,4 +1,4 @@ -from typing import List, Tuple, Optional, Dict, Any +from typing import List, Tuple, Optional, Dict, Any, Union from itertools import combinations, repeat, compress import numpy as np @@ -52,6 +52,47 @@ def __init__( self._fault_planes = fault_planes self._faults = self._calculate_faults(fault_tolerance) + def get_connection_midpoints(self, i: Optional[int] = None) -> np.ndarray: + """ + Returns a numpy array with the midpoint of each connection in the network or, + in case the optional i parameter is specified, the midpoint + of the i'th connection. + + Args: + i: specific zero-offset element to calculate the midpoint for + + Returns: + (Nx3) np.ndarray with connection midpoint coordinates. + + """ + selector: Union[slice, int] = slice(len(self._df_entity_connections.index)) + + if i is not None: + if i > len(self._df_entity_connections.index) or i < 0: + raise ValueError( + f"Optional parameter i is '{i}' but should be between 0 and " + f"{len(self._df_entity_connections.index)}." + ) + if not isinstance(i, int): + raise TypeError( + f"Optional parameter i is of type '{type(i).__name__}' " + "but should be an integer." + ) + selector = i + + coordinates_start = self._df_entity_connections[ + ["xstart", "ystart", "zstart"] + ].values[selector] + coordinates_end = self._df_entity_connections[ + [ + "xend", + "yend", + "zend", + ] + ].values[selector] + + return (coordinates_start + coordinates_end) / 2 + @property def aquifers_xyz(self) -> List[Coordinate]: """ @@ -332,7 +373,7 @@ def _calculate_faults( row["xend"], row["yend"], row["zend"], - *triangle + *triangle, ) if distance: @@ -434,3 +475,28 @@ def df_entity_connections(self) -> pd.DataFrame: def connection_at_nodes(self) -> List[List[int]]: """List is lists with tube indices belonging to a node""" return self._calculate_connections_at_nodes() + + @property + def cell_midpoints(self) -> np.ndarray: + """ + Returns a numpy array with the midpoint of each cell in the network + Returns: + (Nx3) np.ndarray with connection midpoint coordinates. + """ + x_mid = ( + self._grid[["x0", "x1", "x2", "x3", "x4", "x5", "x6", "x7"]] + .mean(axis=1) + .values + ) + y_mid = ( + self._grid[["y0", "y1", "y2", "y3", "y4", "y5", "y6", "y7"]] + .mean(axis=1) + .values + ) + z_mid = ( + self._grid[["z0", "z1", "z2", "z3", "z4", "z5", "z6", "z7"]] + .mean(axis=1) + .values + ) + + return (x_mid, y_mid, z_mid) diff --git a/src/flownet/parameters/_equilibration.py b/src/flownet/parameters/_equilibration.py index b8f942334..18be85d96 100644 --- a/src/flownet/parameters/_equilibration.py +++ b/src/flownet/parameters/_equilibration.py @@ -33,7 +33,7 @@ class Equilibration(Parameter): network: FlowNet network instance. ti2ci: A dataframe with index equal to tube model index, and one column which equals cell indices. eqlnum: A dataframe defining the EQLNUM for each flow tube. - datum_depth: Depth of the datum (m). + datum_depth: Depth of the datum(s) (m). rsvd: Pandas dataframe with a single rsvd table. """ @@ -45,10 +45,10 @@ def __init__( network: NetworkModel, ti2ci: pd.DataFrame, eqlnum: pd.DataFrame, - datum_depth: Optional[float] = None, + datum_depth: Optional[List[float]] = None, rsvd: Optional[pd.DataFrame] = None, ): - self._datum_depth: Union[float, None] = datum_depth + self._datum_depth: Union[List[float], None] = datum_depth self._rsvd: Union[pathlib.Path, None] = rsvd self._network: NetworkModel = network diff --git a/src/flownet/templates/EQUIL.jinja2 b/src/flownet/templates/EQUIL.jinja2 index d896970d6..05f69ec47 100644 --- a/src/flownet/templates/EQUIL.jinja2 +++ b/src/flownet/templates/EQUIL.jinja2 @@ -6,7 +6,7 @@ EQUIL {% for i in range(nr_eqlnum) -%} {% if parameters[i].get("goc_depth", False) -%} {{ parameters[i].get("goc_depth", "1*") }} {{ parameters[i].get("datum_pressure") }} {{ parameters[i].get("owc_depth") }} 1* {{ parameters[i].get("goc_depth", "1*") }} 1* 1* 1* / -{% else -%} -{{ datum_depth }} {{ parameters[i].get("datum_pressure") }} {{ parameters[i].get("owc_depth") }} 1* {{ parameters[i].get("goc_depth", "1*") }} 1* 1* 1* / +{% else %} +{{ datum_depth[i] }} {{ parameters[i].get("datum_pressure") }} {{ parameters[i].get("owc_depth") }} 1* {{ parameters[i].get("goc_depth", "1*") }} 1* 1* 1* / {%- endif -%} {% endfor -%}