From 79f6bf163eb2910a5bae47baab88f7ee45e9c83d Mon Sep 17 00:00:00 2001 From: "Oddvar Lia (ST MSU GEO)" Date: Wed, 11 Dec 2024 10:46:01 +0100 Subject: [PATCH] Add function to combine multiple petro realizations using facies realization Add function to import and export field parameters to and from ERTBOX grid. Added documentation for update_petro_real and the export and import functions. Add test for update_petro_real.py --- docs/export_and_import_fields.rst | 110 +++++ docs/index.rst | 2 + docs/update_petro_real.rst | 157 ++++++ src/fmu/tools/rms/update_petro_real.py | 379 ++++++++++++++ tests/rms/test_update_petro_real.py | 655 +++++++++++++++++++++++++ 5 files changed, 1303 insertions(+) create mode 100644 docs/export_and_import_fields.rst create mode 100644 docs/update_petro_real.rst create mode 100644 src/fmu/tools/rms/update_petro_real.py create mode 100644 tests/rms/test_update_petro_real.py diff --git a/docs/export_and_import_fields.rst b/docs/export_and_import_fields.rst new file mode 100644 index 00000000..0a118470 --- /dev/null +++ b/docs/export_and_import_fields.rst @@ -0,0 +1,110 @@ +rms.export_and_import_field_parameters +======================================= + +The export and import scripts desctibed below supports a workflow where +both facies and petrophysical properties are updated as field parameters in ERT. +The implemented functions in the module *fmu.tools.rms.update_petro_real* are:: + + export_initial_field_parameters + import_updated_field_parameters + +Both the export and import function the same set of petrophysical variables in the +same workflow. Therefore, they both read this information from a configuration file +in YAML format that is also shared with the function *generate_petro_jobs*. + +Export petrophysical parameters as field parameters to be used in ERT FIELD keyword +------------------------------------------------------------------------------------- + +Example of use of the export function + +.. code-block:: python + + from fmu.tools.rms.update_petro_real import export_initial_field_parameters + from fmu.tools.rms.generate_petro_jobs_for_field_update import read_specification_file + + DEBUG_PRINT=True + CONFIG_FILE_NAME_VALYSAR = "../input/config/field_update/generate_petro_jobs_valysar.yml" + CONFIG_FILE_NAME_THERYS = "../input/config/field_update/generate_petro_jobs_therys.yml" + CONFIG_FILE_NAME_VOLON = "../input/config/field_update/generate_petro_jobs_volon.yml" + ERTBOX_GRID = "ERTBOX" + + + def export_fields(): + spec_dict = read_specification_file(CONFIG_FILE_NAME_VALYSAR) + used_petro_dict = spec_dict["used_petro_var"] + export_initial_field_parameters( + project, + used_petro_dict, + grid_model_name=ERTBOX_GRID, + zone_name_for_single_zone_grid="Valysar", + debug_print= DEBUG_PRINT) + + spec_dict = read_specification_file(CONFIG_FILE_NAME_THERYS) + used_petro_dict = spec_dict["used_petro_var"] + export_initial_field_parameters( + project, + used_petro_dict, + grid_model_name=ERTBOX_GRID, + zone_name_for_single_zone_grid="Therys", + debug_print= DEBUG_PRINT) + + spec_dict = read_specification_file(CONFIG_FILE_NAME_VOLON) + used_petro_dict = spec_dict["used_petro_var"] + export_initial_field_parameters( + project, + used_petro_dict, + grid_model_name=ERTBOX_GRID, + zone_name_for_single_zone_grid="Volon", + debug_print= DEBUG_PRINT) + + if __name__ == "__main__": + export_fields() + + +Import petrophysical parameters used as field parameters to be used in ERT FIELD keyword +------------------------------------------------------------------------------------------ + +Example of use of the import function + +.. code-block:: python + + from fmu.tools.rms.update_petro_real import import_updated_field_parameters + from fmu.tools.rms.generate_petro_jobs_for_field_update import read_specification_file + + + DEBUG_PRINT = False + CONFIG_FILE_NAME_VALYSAR = "../input/config/field_update/generate_petro_jobs_valysar.yml" + CONFIG_FILE_NAME_THERYS = "../input/config/field_update/generate_petro_jobs_therys.yml" + CONFIG_FILE_NAME_VOLON = "../input/config/field_update/generate_petro_jobs_volon.yml" + ERTBOX_GRID = "ERTBOX" + + def import_fields(): + spec_dict = read_specification_file(CONFIG_FILE_NAME_VALYSAR) + used_petro_dict = spec_dict["used_petro_var"] + import_updated_field_parameters( + project, + used_petro_dict, + grid_model_name=ERTBOX_GRID, + zone_name_for_single_zone_grid="Valysar", + debug_print=DEBUG_PRINT) + + spec_dict = read_specification_file(CONFIG_FILE_NAME_THERYS) + used_petro_dict = spec_dict["used_petro_var"] + import_updated_field_parameters( + project, + used_petro_dict, + grid_model_name=ERTBOX_GRID, + zone_name_for_single_zone_grid="Therys", + debug_print=DEBUG_PRINT) + + spec_dict = read_specification_file(CONFIG_FILE_NAME_VOLON) + used_petro_dict = spec_dict["used_petro_var"] + import_updated_field_parameters( + project, + used_petro_dict, + grid_model_name=ERTBOX_GRID, + zone_name_for_single_zone_grid="Volon", + debug_print=DEBUG_PRINT) + + if __name__ == "__main__": + import_fields() diff --git a/docs/index.rst b/docs/index.rst index 6caa0ff5..5551ddaa 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -18,6 +18,8 @@ Contents: rename_rms_scripts rms_import_localmodule generate_bw_per_facies + update_petro_real + export_and_import_fields utilities rmsvolumetrics2csv ensembles diff --git a/docs/update_petro_real.rst b/docs/update_petro_real.rst new file mode 100644 index 00000000..b70d1a01 --- /dev/null +++ b/docs/update_petro_real.rst @@ -0,0 +1,157 @@ +rms.update_petro_real +======================= + +When running FMU project where field parameters for both facies and +petrophysical properties are updated in ERT simultaneously, +some adjustments are needed in the RMS project to support this type +of workflow. It is necessary to have one petrosim job per facies. + +The module *update_petro_real* will combine values from realization +of a petrophysical property for individual facies into one common +petrophysical realization by using the facies realization as filter. + +Example describing the workflow +-------------------------------- + +The zone Valysar has four different facies:: + + Floodplain + Channel + Crevasse + Coal + +A petrophysical realization is made for PHIT and KLOGH for each of the three +first facies with property names:: + + Floodplain_PHIT + Floodplain_KLOGH + Channel_PHIT + Channel_KLOGH + Crevasse_PHIT + Crevasse_KLOGH + +Neither PHIT nor KLOGH are updated as field parameters for facies +Coal since they have constant (not spatially varying) values and +may even have 0 specified uncertainty. +The variables VSH and VPHYL are chosen not to be used as field parameters +in ERT in assisted histroy matching and therefore not included here, +but still used in the RMS workflow with the prior values. + +The original realizations created for PHIT and KLOGH are conditioned +to facies. In a workflow where we want to use both facies by applying +the Adaptive PluriGaussian Simulation method (APS) and PHIT and KLOG +for some selected facies as field parameters to be updated by ERT, +we will need the separate PHIT and KLOGH realizations for the selected +facies. + +The steps will be: + + * Use original petrophysical job to create initial version of a realization of + PHIT, KLOGH, VSH and VPHYL. + * Use separate petrophysical jobs to create separate version of realizations of + Floodplain_PHIT, Floodplain_KLOGH and so on. + * Use the script *update_petro_real* to copy the values for PHIT and KLOGH from + the individual realizations per facies into the initial version of PHIT and KLOG + by overwriting the values in the original version. In this process, + the facies realization is used as a filter to select which grid cell values + to get from the individual PHIT and KLOGH parameters that belongs to the + various facies. + +When running ERT iteration 0 which creates the initial ensemble, this looks a bit unnecessary, +since the initial PHIT and KLOGH already has taken facies into account. But, when running +ERT iteration larger than 0 (after ERT has updated the fields Floodplain_PHIT, +Floodplain_KLOGH and so on), then updated versions of the field parameters are imported +and updated version of facies realization is used to update the PHIT and KLOGH parameters +from the imported field parameters (Floodplain_PHIT, Floodplain_KLOGH ...). + +Also for iteration 0, to ensure consistency the same procedure +(copy from the field parameters Floodplain_PHIT, Floodplain_KLOG,..) is applied since this ensure +consistent handling of the realizations. + +Example of RMS python job running this module to combine the field parameters +------------------------------------------------------------------------------- + +The python job below will call the *update_petro_real* function that does the +combination of the field parameters into one petrophysical realization. +This job depends on a small config file ( the same config file that is used by the +function *generate_petro_jobs* from fmu.tools.rms module). Here the python job gets +the grid name and which of the petrophysical properties for which facies +to be updated as field parameters in ERT. The script also depends on the facies table +(which facies code and name belongs together). +This can be fetched from the *global_variables.yml* file as shown in the example below. + +.. code-block:: python + + from fmu.tools.rms.update_petro_real import update_petro_real + from fmu.config import utilities as ut + + DEBUG_PRINT = False + + # Choose either to use config file or input from global variable file + USE_CONFIG_FILES = True + + CFG = ut.yaml_load("../../fmuconfig/output/global_variables.yml")["global"] + FACIES_ZONE = CFG["FACIES_ZONE"] + + # Used if alternative 1 with config file is used + CONFIG_FILES = { + "Valysar": "../input/config/field_update/generate_petro_jobs_valysar.yml", + "Therys": "../input/config/field_update/generate_petro_jobs_therys.yml", + "Volon": "../input/config/field_update/generate_petro_jobs_volon.yml", + } + + # Used if alternative 2 without config file is used + # In this case the global master config must be updated + # to define the dictionary with used petro fields per zone per facies + USED_PETRO_FIELDS = CFG["USED_PETRO_FIELDS"] + GRID_NAMES = { + "Valysar": "Geogrid_Valysar", + "Therys": "Geogrid_Therys", + "Volon": "Geogrid_Volon", + } + FACIES_REAL_NAMES = { + "Valysar": "FACIES", + "Therys": "FACIES", + "Volon": "FACIES", + } + + # For multi zone grids also zone_param_name and zone_code_names must be specified + # Example: + # ZONE_PARAM_NAME = "Zone" + # ZONE_CODE_NAMES = { + # 1: "Valysar, + # 2: "Therys", + # 3: "Volon", + # } + + + if __name__ == "__main__": + # Drogon uses 3 different grids and single zone grids + if USE_CONFIG_FILES: + # Alternative 1 using config file common with generate_petro_jobs + for zone_name in ["Valysar", "Therys", "Volon"]: + facies_code_names = FACIES_ZONE[zone_name] + config_file = CONFIG_FILES[zone_name] + update_petro_real( + project, + facies_code_names, + zone_name_for_single_zone_grid=zone_name, + config_file=config_file, + debug_print=DEBUG_PRINT) + else: + # Alternative 2 specify input using global_variables file + for zone_name in ["Valysar", "Therys", "Volon"]: + facies_code_names = FACIES_ZONE[zone_name] + grid_name = GRID_NAMES[zone_name] + facies_real_name = FACIES_REAL_NAMES[zone_name] + used_petro_per_facies=USED_PETRO_FIELDS + update_petro_real( + project, + facies_code_names, + grid_name=grid_name, + facies_real_name=facies_real_name, + used_petro_dict=used_petro_per_facies, + zone_name_for_single_zone_grid=zone_name, + debug_print=DEBUG_PRINT) + + diff --git a/src/fmu/tools/rms/update_petro_real.py b/src/fmu/tools/rms/update_petro_real.py new file mode 100644 index 00000000..e4544123 --- /dev/null +++ b/src/fmu/tools/rms/update_petro_real.py @@ -0,0 +1,379 @@ +""" +Merge petrophysical realizations created individually per facies +into one realization using facies realization as filter +""" + +from pathlib import Path +from typing import Dict, List + +import xtgeo + +from fmu.tools.rms.generate_petro_jobs_for_field_update import ( + get_original_job_settings, + read_specification_file, +) + + +def update_petro_real( + project, + facies_code_names: Dict[int, str], + config_file: str = "", + grid_name: str = "", + facies_real_name: str = "", + used_petro_dict: Dict = {}, + zone_name_for_single_zone_grid: str = "", + zone_code_names: Dict[int, str] = {}, + zone_param_name: str = "", + debug_print: bool = False, + ignore_missing_parameters: bool = False, +) -> None: + """Combine multiple petrophysical realizations (one per facies) into one parameter + using facies realization as filter. + + Description + + This function will read petrophysical realization for multiple facies + (where all grid cells have the same facies) and use the facies + realization to combine them into one realization conditioned to facies. + + Input + + Choose either to use a config file of the same type as for the function + generate_petro_real in fmu.tools.rms or choose to specify the + dictionary defining which petro variables to use as field parameters for + each zone and facies. In the last case also specify grid model + and facies realization name. + + If multi zone grid is used, also specify dictionary defining + zone name and zone code and zone parameter name. + + Output + + Updated version of petrophysical realizations + for the specified petrophysical variables. + + """ + if config_file: + spec_dict = read_specification_file(config_file) + used_petro_dict = spec_dict["used_petro_var"] + grid_name = spec_dict["grid_name"] + original_job_name = spec_dict["original_job_name"] + + # Get facies param name from the job settings + owner_string_list = ["Grid models", grid_name, "Grid"] + job_type = "Petrophysical Modeling" + petro_job_param = get_original_job_settings( + owner_string_list, job_type, original_job_name + ) + # Get facies realization name from the original petrophysics job + facies_real_name = petro_job_param["InputFaciesProperty"][2] + else: + # Check that necessary parameters are specified + if not grid_name: + raise ValueError( + "Need to specify grid name when config file is not specified." + ) + if not used_petro_dict: + raise ValueError( + "Need to specify the dict used_petro_dict when " + "config file is not specified." + ) + if not facies_real_name: + raise ValueError( + "Need to specify the facies realization name when " + "config file is not specified." + ) + grid = xtgeo.grid_from_roxar(project, grid_name) + subgrids = grid.subgrids + + if subgrids: + # Multi zone grid is found + if not zone_code_names: + raise ValueError( + "Need to specify the 'zone_code_names' when using multi zone grids." + ) + if not zone_param_name: + raise ValueError( + "Need to specify the 'zone_param_name' when using multi zone grids." + ) + if zone_name_for_single_zone_grid: + raise ValueError( + f"For multi zone grids {grid_name} the variable " + "'zone_name_for_single_zone_grid' should not be used" + ) + else: + if not zone_name_for_single_zone_grid: + raise ValueError( + "Need to specify 'zone_name_for_single_zone_grid' " + f"since input {grid_name} is a single zone grid." + ) + + combine_petro_real_from_multiple_facies( + project, + grid_name, + facies_real_name, + used_petro_dict, + facies_code_names, + zone_name_for_single_zone_grid=zone_name_for_single_zone_grid, + zone_param_name=zone_param_name, + zone_code_names=zone_code_names, + debug_print=debug_print, + ignore_missing_parameters=ignore_missing_parameters, + ) + + +def import_updated_field_parameters( + project, + used_petro_dict: Dict, + grid_model_name: str = "ERTBOX", + zone_name_for_single_zone_grid: str = "", + import_path: str = ".", + debug_print: bool = False, +) -> None: + """Import ROFF files with field parameters updated by ERT. + + Description + + This function will import ROFF format files generated by ERT when using + the FIELD keyword in ERT to update petrophysical field parameters. + The naming convention is files with the name of the form: + zonename_faciesname_petrovarname with suffix ".roff" + The files are assumed to be located at the top directory level + where updated fields are written by ERT. + + Input + + A dictionary specifying which petrophysical variables to use as field parameters + for each facies for each zone. + The grid model name to import the field parameters into (ERTBOX grid). + For singe zone grids, also a name of the zone for the single zone + must be specified. + + The result will be new petrophysical parameters in ERTBOX grid. + + """ + + for zone_name, petro_per_facies_dict in used_petro_dict.items(): + if len(used_petro_dict) == 1: + # Single zone grid, use specified zone name + zone_name = zone_name_for_single_zone_grid + if debug_print: + print(f"Zone name: {zone_name}") + for fname, petro_list in petro_per_facies_dict.items(): + if debug_print: + print(f"Facies name: {fname}") + for petro_name in petro_list: + if debug_print: + print(f"Petro variable: {petro_name}") + property_name = zone_name + "_" + fname + "_" + petro_name + file_name = Path(import_path) / Path(property_name + ".roff") + print(f"Import file: {file_name} into {grid_model_name}") + xtgeo_prop = xtgeo.gridproperty_from_file( + file_name, fformat="roff", name=property_name + ) + xtgeo_prop.to_roxar(project, grid_model_name, property_name) + + +def export_initial_field_parameters( + project, + used_petro_dict: Dict, + grid_model_name: str = "ERTBOX", + zone_name_for_single_zone_grid: str = "", + export_path: str = ".", + debug_print: bool = False, +) -> None: + """Export ROFF files with field parameters simulated by RMS to files to be + read by ERT and used as field parameters. + + Description + + This function will export ROFF format files generated by RMS in workflows + where field parameters are updated by ERT. + The parameter names will be in the format: zonename_faciesname_petroname + and the file names will have extension ".roff". + The files are assumed to be located in the directory: rms/output/aps + together with field parameters used by the facies modelling method APS. + + + Input + + A dictionary specifying which petrophysical variables to use as field parameters + for each facies for each zone. + The grid model name to export the field parameters from (ERTBOX grid). + For singe zone grids, also a name of the zone for the single zone must + be specified. + + The result will be files with petrophysical parameters per zone per facies with + size equal to the ERTBOX grid. + + """ + + for zone_name, petro_per_facies_dict in used_petro_dict.items(): + if len(used_petro_dict) == 1: + # Single zone grid, use specified zone name + zone_name = zone_name_for_single_zone_grid + if debug_print: + print(f"Zone name: {zone_name}") + for fname, petro_list in petro_per_facies_dict.items(): + if debug_print: + print(f"Facies name: {fname}") + for petro_name in petro_list: + if debug_print: + print(f"Petro variable: {petro_name}") + property_name = zone_name + "_" + fname + "_" + petro_name + file_name = Path(export_path) / Path(property_name + ".roff") + print(f"Export file: {file_name} into {grid_model_name}") + xtgeo_prop = xtgeo.gridproperty_from_roxar( + project, grid_model_name, property_name + ) + xtgeo_prop.to_file(file_name, fformat="roff", name=property_name) + + +def combine_petro_real_from_multiple_facies( + project, + grid_name: str, + facies_real_name: str, + used_petro_dict: Dict[str, Dict[str, List[str]]], + facies_code_names: Dict[int, str], + zone_name_for_single_zone_grid: str = "", + zone_param_name: str = "", + zone_code_names: Dict[int, str] = {}, + debug_print: bool = False, + ignore_missing_parameters: bool = False, +) -> None: + single_zone_grid = False + if len(zone_name_for_single_zone_grid) > 0: + single_zone_grid = True + + # Find all defined 3D grid parameters using rmsapi + properties = project.grid_models[grid_name].properties + property_names = [prop.name for prop in properties] + # print(f"Gridname: {grid_name} Properties: {property_names}") + + # Find all petro var names to use in any zone + petro_var_list = get_petro_var(used_petro_dict) + + # Get facies realization + prop_facies = xtgeo.gridproperty_from_roxar( + project, grid_name, facies_real_name, faciescodes=True + ) + prop_facies_values = prop_facies.values + + # Get zone realization for multi zone grid + if not single_zone_grid: + prop_zone = xtgeo.gridproperty_from_roxar(project, grid_name, zone_param_name) + prop_zone_values = prop_zone.values + + err_msg = [] + for zone_name, petro_per_facies_dict in used_petro_dict.items(): + if single_zone_grid: + # This is a single zone grid + if len(used_petro_dict) > 1 and zone_name_for_single_zone_grid != zone_name: + # Skip all but the one with correct zone name + continue + # Use the specified zone name + zone_name = zone_name_for_single_zone_grid + else: + if zone_code_names: + zone_code = code_per_name(zone_code_names, zone_name) + for pname in petro_var_list: + prop_petro_original = xtgeo.gridproperty_from_roxar( + project, grid_name, pname + ) + prop_petro_original_values = prop_petro_original.values + is_updated = False + for fname in petro_per_facies_dict: + petro_list_for_this_facies = petro_per_facies_dict[fname] + if pname in petro_list_for_this_facies: + petro_name_per_facies = f"{fname}_{pname}" + + # Get petro realization for this facies and this petro variable + if petro_name_per_facies not in property_names: + err_msg.append( + "Skip non-existing petro realization: " + f"{petro_name_per_facies}" + ) + continue + if ( + project.grid_models[grid_name] + .properties[petro_name_per_facies] + .is_empty() + ): + err_msg.append( + f"Skip empty petro realization: {petro_name_per_facies}" + ) + continue + + prop_petro = xtgeo.gridproperty_from_roxar( + project, grid_name, petro_name_per_facies + ) + prop_petro_values = prop_petro.values + + facies_code = code_per_name(facies_code_names, fname) + + if not single_zone_grid: + # Multi zone grid + if debug_print: + print( + f"Update values for {pname} " + f"in existing parameter for facies {fname} " + f"for zone {zone_name}" + ) + cells_selected = (prop_facies_values == facies_code) & ( + prop_zone_values == zone_code + ) + else: + if debug_print: + print( + f"Update values for {pname} " + f"in existing parameter for facies {fname}" + ) + cells_selected = prop_facies_values == facies_code + is_updated = True + prop_petro_original_values[cells_selected] = prop_petro_values[ + cells_selected + ] + prop_petro_original.values = prop_petro_original_values + if is_updated: + if not single_zone_grid: + # Multi zone grid + print( + f"Write updated petro param {pname} " + f"for zone {zone_name} to grid model {grid_name}" + ) + else: + print( + f"Write updated petro param {pname} for grid model {grid_name}" + ) + prop_petro_original.to_roxar(project, grid_name, pname) + if not ignore_missing_parameters and len(err_msg) > 0: + print( + f"Missing or empty petrophysical 3D parameters for grid model: {grid_name}:" + ) + for msg in err_msg: + print(f" {msg}") + raise ValueError("Missing or empty petrophysical parameters.") + print(f"Finished updating properties for grid model: {grid_name}") + print(" ") + if len(err_msg) > 0: + print( + "Warning: Some petrophysical parameters were not updated. Is that correct?" + ) + + +def code_per_name(code_name_dict: Dict[int, str], input_name: str) -> int: + # Since name is (must be) unique, get it if found or return -1 if not found + for code, name in code_name_dict.items(): + if input_name == name: + return code + return -1 + + +def get_petro_var(used_petro_dict: Dict[str, Dict[str, List[str]]]) -> List[str]: + petro_var_list = [] + for _, petro_var_per_facies_dict in used_petro_dict.items(): + for _, petro_list in petro_var_per_facies_dict.items(): + for petro_name in petro_list: + if petro_name not in petro_var_list: + petro_var_list.append(petro_name) + return petro_var_list diff --git a/tests/rms/test_update_petro_real.py b/tests/rms/test_update_petro_real.py new file mode 100644 index 00000000..3dff8567 --- /dev/null +++ b/tests/rms/test_update_petro_real.py @@ -0,0 +1,655 @@ +"""Run tests in RMS using rmsapi to test update_petro_real.py + +Creates a tmp RMS project in given version. + + +This requires a RMSAPI license, and to be ran in a "roxenvbash" environment + +""" + +import contextlib +import shutil +from os.path import isdir +from pathlib import Path + +import numpy as np +import pytest +import xtgeo # type: ignore + +with contextlib.suppress(ImportError): + import rmsapi # type: ignore + + +from fmu.tools.rms.update_petro_real import ( + export_initial_field_parameters, + import_updated_field_parameters, + update_petro_real, +) + +# ====================================================================================== +# settings to create RMS project! +DEBUG_PRINT = False +REMOVE_RMS_PROJECT_AFTER_TEST = True + +TMPD = Path("TMP") +TMPD.mkdir(parents=True, exist_ok=True) + + +PROJNAME = "tmp_project_update_petro_real.rmsxxx" +PRJ = str(TMPD / PROJNAME) +RESULTDIR = TMPD / "update_petro" +RESULTDIR.mkdir(parents=True, exist_ok=True) + +EXPORT_PATH = RESULTDIR +IMPORT_PATH = RESULTDIR + +FACIES_CODE_NAMES = { + 1: "F1", + 2: "F2", + 3: "F3", + 4: "F4", +} + +ZONE_CODE_NAMES = { + 1: "ZoneA", + 2: "ZoneB", +} + + +USED_PETRO_DICT = { + "ZoneA": { + "F1": ["P1"], + "F2": ["P1", "P2"], + }, + "ZoneB": { + "F3": ["P2"], + "F2": ["P1", "P2"], + }, +} +# This grid is always a single zone grid +GRID_MODEL_ERTBOX = "ERTBOX" + +# This grid is a single zone geogrid +GRID_MODEL_GEO_SINGLE = "SingleZoneGrid" + +# This grid is a multi zone geogrid +GRID_MODEL_GEO_MULTI = "MultiZoneGrid" + +FACIES_REAL_NAME = "Facies" +ZONE_PARAM_NAME = "Zones" + +NX = 5 +NY = 6 +NZ_ERTBOX = 5 +NZ = 10 + +SUBGRID_DICT = { + "ZoneA": int(NZ / 2), + "ZoneB": NZ - int(NZ / 2), +} + +ZONE_FACIES_PETRO_VALUES = { + "ZoneA": { + "F1": { + "P1": 10.0, + "P2": 100.0, + }, + "F2": { + "P1": 20.0, + "P2": 150.0, + }, + "F3": { + "P1": 14, + "P2": 130.0, + }, + "F4": { + "P1": 28, + "P2": 165.0, + }, + }, + "ZoneB": { + "F1": { + "P1": 16.0, + "P2": 110.0, + }, + "F2": { + "P1": 22.0, + "P2": 160.0, + }, + "F3": { + "P1": 12.0, + "P2": 120.0, + }, + "F4": { + "P1": 25.0, + "P2": 175.0, + }, + }, +} + + +def make_facies_param_multizone(dimensions: tuple[int, int, int]): + values = np.zeros(dimensions, dtype=np.uint8) + sum_layer = 0 + nfacies = len(FACIES_CODE_NAMES) + for zone_name, nlayers in SUBGRID_DICT.items(): + start_layer = sum_layer + end_layer = start_layer + nlayers + petro_per_facies_dict = USED_PETRO_DICT[zone_name] + facies_names = list(petro_per_facies_dict.keys()) + nfacies = len(facies_names) + for layer in range(start_layer, end_layer): + indx = layer % nfacies + fname = facies_names[indx] + facies_code = get_facies_code(fname) + values[:, :, layer] = facies_code + sum_layer += nlayers + return values + + +def make_facies_param_singlezone(dimensions: tuple[int, int, int]): + values = np.zeros(dimensions, dtype=np.uint8) + for k in range(dimensions[2]): + facies_for_layer = (k % 2) + 1 + values[:, :, k] = facies_for_layer + return values + + +def make_zone_param(dimensions: tuple[int, int, int]): + values = np.zeros(dimensions, dtype=np.uint8) + sum_layer = 0 + for zone_name, nlayers in SUBGRID_DICT.items(): + start_layer = sum_layer + end_layer = start_layer + nlayers + values[:, :, start_layer:end_layer] = get_zone_code(zone_name) + sum_layer += nlayers + return values + + +def make_petro_param(dimensions: tuple[int, int, int], value: float): + values = np.zeros(dimensions, dtype=np.float32) + values[:, :, :] = value + return values + + +def make_reference_petro_param( + dimensions: tuple[int, int, int], petro_name: str, facies_values, zone_values +): + (nx, ny, nz) = dimensions + values = np.zeros(dimensions, dtype=np.float32) + # Assign only values for the grid cells containing wanted facies and zone code + # All other grid cells are kept to 0. They will not be updated anyway and it + # simplifies the test if they are 0. + for k in range(nz): + for j in range(ny): + for i in range(nx): + facies_code = facies_values[i, j, k] + facies_name = FACIES_CODE_NAMES[facies_code] + zone_code = zone_values[i, j, k] + zone_name = ZONE_CODE_NAMES[zone_code] + if zone_name in USED_PETRO_DICT: + petro_per_facies_dict = USED_PETRO_DICT[zone_name] + if facies_name in petro_per_facies_dict: + petro_list = petro_per_facies_dict[facies_name] + if petro_name in petro_list: + values[i, j, k] = ZONE_FACIES_PETRO_VALUES[zone_name][ + facies_name + ][petro_name] + return values + + +def get_zone_code(zone_name: str): + for zone_code, name in ZONE_CODE_NAMES.items(): + if name == zone_name: + return zone_code + raise ValueError(f"Can not find zone {zone_name} in 'ZONE_CODE_NAMES'") + + +def get_facies_code(facies_name: str): + for facies_code, name in FACIES_CODE_NAMES.items(): + if name == facies_name: + return facies_code + raise ValueError(f"Can not find facies {facies_name} in 'FACIES_CODE_NAMES'") + + +def make_petro_param_per_zone( + dimensions: tuple[int, int, int], + input_values, + zone_name: str, + zone_values, + constant_value: float, +): + # Will modify input_values and return updated values + (nx, ny, nz) = dimensions + if input_values is None: + values = np.zeros(dimensions, dtype=np.float32) + else: + values = input_values + zone_code = get_zone_code(zone_name) + selected = zone_values == zone_code + values[selected] = constant_value + return values + + +def get_petro_variable_names(): + petro_names = [] + for _, petro_per_facies_dict in USED_PETRO_DICT.items(): + for _, petro_list in petro_per_facies_dict.items(): + for petro_name in petro_list: + if petro_name not in petro_names: + petro_names.append(petro_name) + return petro_names + + +def create_project(): + """Create a tmp RMS project for testing, populate with basic data.""" + prj1 = str(PRJ) + + print("\n******** Setup RMS project!\n") + if isdir(prj1): + print("Remove existing project! (1)") + shutil.rmtree(prj1) + + project = rmsapi.Project.create() + + rox = xtgeo.RoxUtils(project) + print("rmsapi version is", rox.roxversion) + print("RMS version is", rox.rmsversion(rox.roxversion)) + assert "1." in rox.roxversion + + dimensions = (NX, NY, NZ_ERTBOX) + xtgeo_ertbox_grid = xtgeo.create_box_grid(dimensions, increment=(50.0, 50.0, 1.0)) + xtgeo_ertbox_grid.to_roxar(project, GRID_MODEL_ERTBOX) + + dimensions = (NX, NY, NZ) + xtgeo_grid_single = xtgeo.create_box_grid(dimensions, increment=(50.0, 50.0, 1.0)) + xtgeo_grid_single.to_roxar(project, GRID_MODEL_GEO_SINGLE) + + xtgeo_grid_multizone = xtgeo.create_box_grid( + dimensions, increment=(50.0, 50.0, 1.0) + ) + xtgeo_grid_multizone.set_subgrids(SUBGRID_DICT) + xtgeo_grid_multizone.to_roxar(project, GRID_MODEL_GEO_MULTI) + + values = make_facies_param_multizone((NX, NY, NZ)) + xtgeo_facies_param_multi = xtgeo.GridProperty( + ncol=NX, + nrow=NY, + nlay=NZ, + name=FACIES_REAL_NAME, + discrete=True, + grid=xtgeo_grid_multizone, + codes=FACIES_CODE_NAMES, + roxar_dtype=np.uint8, + values=values, + ) + xtgeo_facies_param_multi.to_roxar(project, GRID_MODEL_GEO_MULTI, FACIES_REAL_NAME) + + values = make_facies_param_singlezone((NX, NY, NZ)) + xtgeo_facies_param_single = xtgeo.GridProperty( + ncol=NX, + nrow=NY, + nlay=NZ, + name=FACIES_REAL_NAME, + discrete=True, + grid=xtgeo_grid_single, + codes=FACIES_CODE_NAMES, + roxar_dtype=np.uint8, + values=values, + ) + xtgeo_facies_param_single.to_roxar(project, GRID_MODEL_GEO_SINGLE, FACIES_REAL_NAME) + + values = make_zone_param((NX, NY, NZ)) + xtgeo_zone_param_multi = xtgeo.GridProperty( + ncol=NX, + nrow=NY, + nlay=NZ, + name=ZONE_PARAM_NAME, + discrete=True, + grid=xtgeo_grid_multizone, + codes=ZONE_CODE_NAMES, + roxar_dtype=np.uint8, + values=values, + ) + xtgeo_zone_param_multi.to_roxar(project, GRID_MODEL_GEO_MULTI, ZONE_PARAM_NAME) + + values = np.ones((NX, NY, NZ), dtype=np.uint8) + xtgeo_zone_param_single = xtgeo.GridProperty( + ncol=NX, + nrow=NY, + nlay=NZ, + name=ZONE_PARAM_NAME, + discrete=True, + grid=xtgeo_grid_single, + codes=ZONE_CODE_NAMES, + roxar_dtype=np.uint8, + values=values, + ) + xtgeo_zone_param_single.to_roxar(project, GRID_MODEL_GEO_SINGLE, ZONE_PARAM_NAME) + + # make petro param per zone and facies for ertbox grid + for zone_name, petro_per_facies_dict in USED_PETRO_DICT.items(): + for facies_name, petro_list in petro_per_facies_dict.items(): + for petro_name in petro_list: + new_petro_name = zone_name + "_" + facies_name + "_" + petro_name + # Use constant value for all grid cells + values = make_petro_param( + (NX, NY, NZ_ERTBOX), + ZONE_FACIES_PETRO_VALUES[zone_name][facies_name][petro_name], + ) + xtgeo_petro_param = xtgeo.GridProperty( + ncol=NX, + nrow=NY, + nlay=NZ_ERTBOX, + name=new_petro_name, + grid=xtgeo_ertbox_grid, + roxar_dtype=np.float32, + values=values, + ) + xtgeo_petro_param.to_roxar(project, GRID_MODEL_ERTBOX, new_petro_name) + + # make inital and reference params for petro variables for geo grids + petro_names = get_petro_variable_names() + for petro_name in petro_names: + petro_name_ref = petro_name + "_ref" + facies_values_multi = xtgeo_facies_param_multi.values + zone_values_multi = xtgeo_zone_param_multi.values + + facies_values_single = xtgeo_facies_param_single.values + zone_values_single = xtgeo_zone_param_single.values + + # The reference petro params will depend on zone, + # facies and be constant for given zone and facies + values = make_reference_petro_param( + (NX, NY, NZ), petro_name, facies_values_multi, zone_values_multi + ) + xtgeo_petro_param = xtgeo.GridProperty( + ncol=NX, + nrow=NY, + nlay=NZ, + name=petro_name_ref, + grid=xtgeo_grid_multizone, + roxar_dtype=np.float32, + values=values, + ) + xtgeo_petro_param.to_roxar(project, GRID_MODEL_GEO_MULTI, petro_name_ref) + + values = make_reference_petro_param( + (NX, NY, NZ), petro_name, facies_values_single, zone_values_single + ) + xtgeo_petro_param = xtgeo.GridProperty( + ncol=NX, + nrow=NY, + nlay=NZ, + name=petro_name_ref, + grid=xtgeo_grid_single, + roxar_dtype=np.float32, + values=values, + ) + xtgeo_petro_param.to_roxar(project, GRID_MODEL_GEO_SINGLE, petro_name_ref) + + # Set initial value of petro params to 0 + # They will be updated and checked with reference after that. + values = np.zeros((NX, NY, NZ), dtype=np.float32) + xtgeo_petro_param_initial = xtgeo.GridProperty( + ncol=NX, + nrow=NY, + nlay=NZ, + name=petro_name_ref, + grid=xtgeo_grid_single, + roxar_dtype=np.float32, + values=values, + ) + xtgeo_petro_param_initial.to_roxar(project, GRID_MODEL_GEO_MULTI, petro_name) + xtgeo_petro_param_initial.to_roxar(project, GRID_MODEL_GEO_SINGLE, petro_name) + + # make petro param per facies for multizone grid. + # They are use as input when updating the petro params + values = None + zone_values_multi = xtgeo_zone_param_multi.values + values_petro_param_dict = {} + for zone_name, petro_per_facies_dict in USED_PETRO_DICT.items(): + for facies_name, petro_list in petro_per_facies_dict.items(): + if facies_name not in values_petro_param_dict: + values_petro_param_dict[facies_name] = {} + for petro_name in petro_list: + if petro_name not in values_petro_param_dict[facies_name]: + # Initialize all values to 0 + values = np.zeros((NX, NY, NZ), dtype=np.float32) + values_petro_param_dict[facies_name][petro_name] = values + else: + values = values_petro_param_dict[facies_name][petro_name] + + # Update values for current zone, facies and petro_param + # with a constant value + values_petro_param_dict[facies_name][petro_name] = ( + make_petro_param_per_zone( + (NX, NY, NZ), + values, + zone_name, + zone_values_multi, + ZONE_FACIES_PETRO_VALUES[zone_name][facies_name][petro_name], + ) + ) + + new_petro_param_list = [] + for zone_name, petro_per_facies_dict in USED_PETRO_DICT.items(): + for facies_name, petro_list in petro_per_facies_dict.items(): + for petro_name in petro_list: + new_petro_name = facies_name + "_" + petro_name + if new_petro_name not in new_petro_param_list: + new_petro_param_list.append(new_petro_name) + values = values_petro_param_dict[facies_name][petro_name] + xtgeo_petro_param = xtgeo.GridProperty( + ncol=NX, + nrow=NY, + nlay=NZ, + name=new_petro_name, + grid=xtgeo_grid_multizone, + roxar_dtype=np.float32, + values=values, + ) + xtgeo_petro_param.to_roxar( + project, GRID_MODEL_GEO_MULTI, new_petro_name + ) + + # make petro param per facies for single zone grid. + # They are use as input when updating the petro params + values = None + zone_values_single = xtgeo_zone_param_single.values + zone_name = ZONE_CODE_NAMES[1] + petro_per_facies_dict = USED_PETRO_DICT[zone_name] + for facies_name, petro_list in petro_per_facies_dict.items(): + for petro_name in petro_list: + new_petro_name = facies_name + "_" + petro_name + # Constant value for all grid cells in each zone. + # The input 'values' will initially be created here + # and updated for each zone + values = make_petro_param_per_zone( + (NX, NY, NZ), + values, + zone_name, + zone_values_single, + ZONE_FACIES_PETRO_VALUES[zone_name][facies_name][petro_name], + ) + xtgeo_petro_param = xtgeo.GridProperty( + ncol=NX, + nrow=NY, + nlay=NZ, + name=new_petro_name, + grid=xtgeo_grid_single, + roxar_dtype=np.float32, + values=values, + ) + xtgeo_petro_param.to_roxar(project, GRID_MODEL_GEO_SINGLE, new_petro_name) + + project.save_as(prj1) + return project + + +PROJECT = create_project() + + +@pytest.mark.skipunlessroxar +@pytest.mark.parametrize( + "used_petro_dict,zone_name_for_single_zone_grid", + [ + ( + { + "default": USED_PETRO_DICT["ZoneA"], + }, + "ZoneA", + ), + ( + { + "default": USED_PETRO_DICT["ZoneB"], + }, + "ZoneB", + ), + ], +) +def test_import_export_updated_field_parameters( + used_petro_dict: dict, + zone_name_for_single_zone_grid: str, +) -> None: + export_initial_field_parameters( + PROJECT, + used_petro_dict, + grid_model_name=GRID_MODEL_ERTBOX, + zone_name_for_single_zone_grid=zone_name_for_single_zone_grid, + export_path=EXPORT_PATH, + debug_print=DEBUG_PRINT, + ) + + # Compare input with reference data + for zone_name, petro_per_facies_dict in used_petro_dict.items(): + if zone_name == zone_name_for_single_zone_grid: + for facies_name, petro_list in petro_per_facies_dict.items(): + for petro_name in petro_list: + name = zone_name + "_" + facies_name + "_" + petro_name + name_from_export = name + "_" + "export" + filename = RESULTDIR / Path(name + ".roff") + xtgeo_property = xtgeo.gridproperty_from_file( + filename, fformat="roff" + ) + xtgeo_property.to_roxar( + PROJECT, GRID_MODEL_ERTBOX, name_from_export + ) + value1 = ( + PROJECT.grid_models[GRID_MODEL_ERTBOX] + .properties[name_from_export] + .get_values() + ) + value2 = ( + PROJECT.grid_models[GRID_MODEL_ERTBOX] + .properties[name] + .get_values() + ) + assert np.allclose(value1, value2) + + import_updated_field_parameters( + PROJECT, + used_petro_dict, + grid_model_name=GRID_MODEL_ERTBOX, + zone_name_for_single_zone_grid=zone_name_for_single_zone_grid, + import_path=IMPORT_PATH, + debug_print=DEBUG_PRINT, + ) + + # Compare input with reference data + for zone_name, petro_per_facies_dict in used_petro_dict.items(): + if zone_name == zone_name_for_single_zone_grid: + for facies_name, petro_list in petro_per_facies_dict.items(): + for petro_name in petro_list: + name = zone_name + "_" + facies_name + "_" + petro_name + name_from_export = name + "_" + "export" + filename = RESULTDIR / Path(name + ".roff") + xtgeo_property = xtgeo.gridproperty_from_file( + filename, fformat="roff" + ) + xtgeo_property.to_roxar( + PROJECT, GRID_MODEL_ERTBOX, name_from_export + ) + values1 = ( + PROJECT.grid_models[GRID_MODEL_ERTBOX] + .properties[name_from_export] + .get_values() + ) + values2 = ( + PROJECT.grid_models[GRID_MODEL_ERTBOX] + .properties[name] + .get_values() + ) + assert np.allclose(values1, values2) + + +@pytest.mark.skipunlessroxar +def test_update_field_parameters_single_zone_grid() -> None: + update_petro_real( + PROJECT, + facies_code_names=FACIES_CODE_NAMES, + grid_name=GRID_MODEL_GEO_SINGLE, + facies_real_name=FACIES_REAL_NAME, + used_petro_dict=USED_PETRO_DICT, + zone_name_for_single_zone_grid=ZONE_CODE_NAMES[1], + debug_print=DEBUG_PRINT, + ) + + # Compare with reference + petro_names = get_petro_variable_names() + for petro_name in petro_names: + petro_name_reference = petro_name + "_ref" + values1 = ( + PROJECT.grid_models[GRID_MODEL_GEO_SINGLE] + .properties[petro_name] + .get_values() + ) + values2 = ( + PROJECT.grid_models[GRID_MODEL_GEO_SINGLE] + .properties[petro_name_reference] + .get_values() + ) + assert np.allclose(values1, values2) + + +@pytest.mark.skipunlessroxar +def test_update_field_parameters_multi_zone_grid() -> None: + update_petro_real( + PROJECT, + facies_code_names=FACIES_CODE_NAMES, + grid_name=GRID_MODEL_GEO_MULTI, + facies_real_name=FACIES_REAL_NAME, + used_petro_dict=USED_PETRO_DICT, + zone_code_names=ZONE_CODE_NAMES, + zone_param_name=ZONE_PARAM_NAME, + debug_print=DEBUG_PRINT, + ) + + # Compare with reference + petro_names = get_petro_variable_names() + for petro_name in petro_names: + petro_name_reference = petro_name + "_ref" + values1 = ( + PROJECT.grid_models[GRID_MODEL_GEO_MULTI] + .properties[petro_name] + .get_values() + ) + values2 = ( + PROJECT.grid_models[GRID_MODEL_GEO_MULTI] + .properties[petro_name_reference] + .get_values() + ) + assert np.allclose(values1, values2) + + +@pytest.mark.skipunlessroxar +def test_close_project(): + PROJECT.save() + if REMOVE_RMS_PROJECT_AFTER_TEST: + print("\n******* Teardown RMS project!\n") + if isdir(PRJ): + print("Remove existing project!") + shutil.rmtree(PRJ) + if isdir(RESULTDIR): + print("Remove temporary files") + shutil.rmtree(RESULTDIR)