From d9dae16ed1c99296540e96c440013be558af4ce5 Mon Sep 17 00:00:00 2001 From: "Oddvar Lia (ST MSU GEO)" Date: Fri, 20 Oct 2023 10:38:42 +0200 Subject: [PATCH] update sim_field.py --- .../example_case/scripts/__init__.py | 0 .../example_case/scripts/common_functions.py | 443 +++++++++++ .../example_case/scripts/init_test_case.py | 283 +++++++ .../example_case/scripts/sim_fields.py | 701 ++---------------- 4 files changed, 778 insertions(+), 649 deletions(-) create mode 100644 tests/jobs/localisation/example_case/scripts/__init__.py create mode 100644 tests/jobs/localisation/example_case/scripts/common_functions.py create mode 100755 tests/jobs/localisation/example_case/scripts/init_test_case.py diff --git a/tests/jobs/localisation/example_case/scripts/__init__.py b/tests/jobs/localisation/example_case/scripts/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/jobs/localisation/example_case/scripts/common_functions.py b/tests/jobs/localisation/example_case/scripts/common_functions.py new file mode 100644 index 000000000..27366aefd --- /dev/null +++ b/tests/jobs/localisation/example_case/scripts/common_functions.py @@ -0,0 +1,443 @@ +""" +Common functions used by the scripts: init_test_case.py and sim_field.py +""" + +import xtgeo +import gaussianfft as sim +import numpy as np + + +# NOTE: xtgeo MUST be imported BEFORE gaussianfft +# The reason is that xtgeo has functions importing roxar API and +# even though this is not used, the code will crash with core dump +# since Roxar API and gaussianfft both use Boost to wrap C++ code +# into python functions but Roxar API and gaussianfft uses two +# slightly different versions of Boost. +# The gaussianfft module uses version 1.76 which is newer than +# version 1.74 from Roxar API (and indirectly xtgeo) +# which may explain why it works importing gaussianfft after +# xtgeo (and Roxar API module roxar). We don't know exactly the reason +# for the core dump with wrong sequence of the import's, but probably +# it is due to some initialization related to the Boost library and +# Boost version 1.74.0 is not compatible with any initialization done by +# version 1.76. +# +# The message related to Boost and RMS (and then indirectly also xtgeo) +# from Aspentech Support is this: +# "The current version of boost is 1.74.0 and this version has been used +# since RMS 12.1, and is still used in RMS V14.0.1 and V14.1. +# Boost version 1.81.0 will be available in version 14.2 or version 15. +# Thomas also write "Last time I checked, boost did not provide any +# compatibility guarantees, so it's not expected to work if you +# mix two different boost versions in the same process +# by loading python modules into RMS that uses other versions." + +# The current combination of gaussianfft and xtgeo (or RMS Roxar API) +# will work if xtgeo is imported first. But this may change later. +# The plan for gaussianfft is to ensure correct sequence of import by +# importing xtgeo (and when running from RMS, also roxar) in +# gaussianfft before calling any functions from the gaussianfft module. +# In this case it should work for the end user regardless of which sequence +# it is imported. The most robust solution would be to generate gaussianfft +# to use exactly the same. + + +# pylint: disable=missing-function-docstring, too-many-locals, invalid-name +# pylint: disable=bare-except, raise-missing-from +# pylint: disable= redefined-outer-name, too-many-nested-blocks + + +def specify_settings(): + """ + grid_size - length, width, thickness of a box containing the field + Same size is used for both fine scale grid with the simulated field + and the coarse scale grid containing upscaled values + of the simulated field. + + field_settings - Define the dimension (number of grid cells) for fine scale grid, + name of output files and specification of model parameters for + simulation of gaussian field with option to use linear trend. + Relative standard deviation specify standard deviation of + gaussian residual field relative to the trends span of value + (max trend value - min trend value) + + response_settings - Specify the coarse grid dimensions, name of file and type + of average operation to calculated upscaled values that + are predictions of observations of the same grid cells. + Which cell indices are observed are specified in + observation settings. + + observation_settings - Specify name of files for generated observations + and also which grid cells from coarse grid is used + as observables. + (Cells that have values that are used as observations) + """ + grid_size = { + "xsize": 7500.0, + "ysize": 12500.0, + "zsize": 50.0, + "use_eclipse_grid_index_origin": True, + } + field_settings = { + "grid_dimension": [150, 250, 1], + "grid_file_name": "GRID.EGRID", + "field_correlation_range": [3000.0, 1000.0, 2.0], + "field_correlation_azimuth": 45.0, + "field_correlation_dip": 0.0, + "field_variogram": "spherical", + "field_trend_params": [1.0, -1.0], + "field_trend_relstd": 0.05, + "field_trend_use": 0, + "field_name": "FIELDPARAM", + "field_updated_file_name": "FieldParam.roff", + "field_initial_file_name": "init_files/FieldParam.roff", + "field_seed_file": "randomseeds.txt", + } + response_settings = { + "grid_dimension": [15, 25, 1], + "upscaled_file_name": "Upscaled.roff", + "grid_file_name": "UpscaleGrid.EGRID", + "response_function": "average", + "gen_data_file_name": "UpscaledField_0.txt", + "all": True, + } + + observation_settings = { + "observation_dir": "observations", + "observation_file": "observations.obs", + "observation_data_dir": "obs_data", + "3D_param_file_name": "init_files/UpscaledObsField.roff", + "rel_error": 0.10, + "min_abs_error": 0.01, + "selected_grid_cells": [ + [5, 10, 1], + ], + # "selected_grid_cells":[ + # [15, 1, 1], + # [15, 25, 1], + # [ 1, 25, 1], + # [ 1, 1, 1], + # [ 3, 3, 1], + # [ 5, 23, 1], + # [11, 4, 1], + # [ 3, 12, 1], + # [13, 18, 1], + # ], + } + return grid_size, field_settings, response_settings, observation_settings + + +def generate_field_and_upscale( + field_settings, grid_size, response_settings, observation_settings, real_number +): + seed_file_name = field_settings["field_seed_file"] + relative_std = field_settings["field_trend_relstd"] + use_trend = field_settings["field_trend_use"] + + start_seed = get_seed(seed_file_name, real_number) + residual_field = simulate_field(grid_size, field_settings, start_seed) + if use_trend == 1: + trend_field = trend(grid_size, field_settings) + field3D = trend_field + relative_std * residual_field + else: + field3D = residual_field + + # Write field parameter for fine scale grid + field_object = export_field(field_settings, field3D) + field_values = field_object.values + + # Calculate upscaled values for selected coarse grid cells + upscaled_values = upscaling( + field_values, + response_settings, + observation_settings, + write_field=True, + iteration=0, + ) + return upscaled_values + + +def get_seed(seed_file_name, r_number): + with open(seed_file_name, "r", encoding="utf8") as file: + lines = file.readlines() + try: + seed_value = int(lines[r_number - 1]) + except: # noqa: E722 + raise IOError("Seed value not found for realization {r_number} ") + return seed_value + + +def upscaling( + field_values, response_settings, observation_settings, write_field=True, iteration=0 +): + response_function_name = response_settings["response_function"] + upscaled_file_name = response_settings["upscaled_file_name"] + NX, NY, NZ = response_settings["grid_dimension"] + calculate_all = response_settings["all"] + + coarse_cell_index_list = observation_settings["selected_grid_cells"] + upscaled_values = np.zeros((NX, NY, NZ), dtype=np.float32, order="F") + upscaled_values[:, :, :] = -999 + + if response_function_name == "average": + upscaled_values = upscale_average( + field_values, coarse_cell_index_list, upscaled_values, use_all=calculate_all + ) + + if iteration == 0: + upscaled_file_name = "init_files/" + upscaled_file_name + + if write_field: + write_upscaled_field(upscaled_values, upscaled_file_name) + return upscaled_values + + +def write_upscaled_field( + upscaled_values, upscaled_file_name, selected_cell_index_list=None +): + nx, ny, nz = upscaled_values.shape + field_name = "Upscaled" + + field_object = xtgeo.grid3d.GridProperty( + ncol=nx, + nrow=ny, + nlay=nz, + values=upscaled_values, + discrete=False, + name=field_name, + ) + + print(f"Write upscaled field file: {upscaled_file_name} ") + field_object.to_file(upscaled_file_name, fformat="roff") + if selected_cell_index_list is not None: + selected_upscaled_values = np.zeros((nx, ny, nz), dtype=np.float32, order="F") + selected_upscaled_values[:, :, :] = -1 + for indices in selected_cell_index_list: + Iindx = indices[0] - 1 + Jindx = indices[1] - 1 + Kindx = indices[2] - 1 + selected_upscaled_values[Iindx, Jindx, Kindx] = upscaled_values[ + Iindx, Jindx, Kindx + ] + + field_name_selected = field_name + "_conditioned_cells" + file_name_selected = "init_files/" + field_name_selected + ".roff" + cond_field_object = xtgeo.grid3d.GridProperty( + ncol=nx, + nrow=ny, + nlay=nz, + values=selected_upscaled_values, + discrete=False, + name=field_name_selected, + ) + print(f"Write conditioned cell values as field: {file_name_selected}") + cond_field_object.to_file(file_name_selected, fformat="roff") + + return field_object + + +def upscale_average( + field_values, coarse_cell_index_list, upscaled_values, use_all=False +): + """ + Input: field_values (numpy 3D) + coarse_cell_index_list (list of tuples (I,J,K)) + Output: upscaled_values (numpy 3D) initialized outside + but filled in specified (I,J,K) cells. + """ + nx, ny, nz = field_values.shape + NX, NY, NZ = upscaled_values.shape + + print(f"Number of fine scale grid cells: (nx,ny,nz): ({nx},{ny},{nz})") + print(f"Number of coarse scale grid cells: (NX,NY,NZ): ({NX},{NY},{NZ}) ") + mx = int(nx / NX) + my = int(ny / NY) + mz = int(nz / NZ) + print( + "Number of fine scale grid cells per coarse grid cell: " + f"(mx,my,mz): ({mx},{my},{mz}) " + ) + if use_all: + print("Calculate upscaled values for all grid cells") + for Kindx in range(NZ): + for Jindx in range(NY): + for Iindx in range(NX): + istart = mx * Iindx + iend = istart + mx + jstart = my * Jindx + jend = jstart + my + kstart = mz * Kindx + kend = kstart + mz + sum_val = 0.0 + for k in range(kstart, kend): + for j in range(jstart, jend): + for i in range(istart, iend): + sum_val += field_values[i, j, k] + upscaled_values[Iindx, Jindx, Kindx] = sum_val / (mx * my * mz) + + else: + print("Calculate upscaled values for selected grid cells") + for indices in coarse_cell_index_list: + Iindx = indices[0] - 1 + Jindx = indices[1] - 1 + Kindx = indices[2] - 1 + istart = mx * Iindx + iend = istart + mx + jstart = my * Jindx + jend = jstart + my + kstart = mz * Kindx + kend = kstart + mz + sum_val = 0.0 + for k in range(kstart, kend): + for j in range(jstart, jend): + for i in range(istart, iend): + sum_val += field_values[i, j, k] + upscaled_values[Iindx, Jindx, Kindx] = sum_val / (mx * my * mz) + + return upscaled_values + + +def write_prediction_gen_data(upscaled_values, observation_settings, response_settings): + cell_indx_list = observation_settings["selected_grid_cells"] + response_file_name = response_settings["gen_data_file_name"] + print(f"Write GEN_DATA file with prediction of observations: {response_file_name}") + with open(response_file_name, "w", encoding="utf8") as file: + # NOTE: The sequence of values must be the same as for the observations + for indices in cell_indx_list: + Iindx = indices[0] - 1 + Jindx = indices[1] - 1 + Kindx = indices[2] - 1 + value = upscaled_values[Iindx, Jindx, Kindx] + print(f"Prediction of obs for {Iindx+1},{Jindx+1},{Kindx+1}: {value}") + file.write(f"{value}\n") + + +def trend(grid_size, field_settings): + """ + Return 3D numpy array with values following a linear trend + scaled to take values between 0 and 1. + """ + nx, ny, nz = field_settings["grid_dimension"] + xsize = grid_size["xsize"] + ysize = grid_size["ysize"] + a, b = field_settings["field_trend_params"] + + x0 = 0.0 + y0 = 0.0 + dx = xsize / nx + dy = ysize / ny + + maxsize = ysize + if xsize > ysize: + maxsize = xsize + + val = np.zeros((nx, ny, nz), dtype=np.float32, order="F") + for i in range(nx): + x = x0 + i * dx + for j in range(ny): + y = y0 + j * dy + for k in range(nz): + val[i, j, k] = a * (x - x0) / maxsize + b * (y - y0) / maxsize + + minval = np.min(val) + maxval = np.max(val) + val_normalized = (val - minval) / (maxval - minval) + return val_normalized + + +def simulate_field(grid_size, field_settings, start_seed): + # pylint: disable=no-member, + variogram_name = field_settings["field_variogram"] + corr_ranges = field_settings["field_correlation_range"] + xrange = corr_ranges[0] + yrange = corr_ranges[1] + zrange = corr_ranges[2] + + azimuth = field_settings["field_correlation_azimuth"] + dip = field_settings["field_correlation_dip"] + + nx, ny, nz = field_settings["grid_dimension"] + xsize = grid_size["xsize"] + ysize = grid_size["ysize"] + zsize = grid_size["zsize"] + + dx = xsize / nx + dy = ysize / ny + dz = zsize / nz + + print(f"Start seed: {start_seed}") + sim.seed(start_seed) + + variogram = sim.variogram( + variogram_name, + xrange, + perp_range=yrange, + depth_range=zrange, + azimuth=azimuth - 90, + dip=dip, + ) + + print(f"Simulate field with size: nx={nx},ny={ny} ") + field1D = sim.simulate(variogram, nx, dx, ny, dy, nz, dz) + field = field1D.reshape((nx, ny, nz), order="F") + return field + + +def export_field(field_settings, field3D): + # Export initial ensemble field + nx, ny, nz = field_settings["grid_dimension"] + field_name = field_settings["field_name"] + field_file_name = field_settings["field_initial_file_name"] + + field_object = xtgeo.grid3d.GridProperty( + ncol=nx, nrow=ny, nlay=nz, values=field3D, discrete=False, name=field_name + ) + + print(f"Write field file: {field_file_name} ") + field_object.to_file(field_file_name, fformat="roff") + return field_object + + +def read_field_from_file(field_settings): + input_file_name = field_settings["field_updated_file_name"] + name = field_settings["field_name"] + field_object = xtgeo.gridproperty_from_file( + input_file_name, fformat="roff", name=name + ) + return field_object + + +def read_obs_field_from_file(observation_settings): + input_file_name = observation_settings["3D_param_file_name"] + obs_field_object = xtgeo.gridproperty_from_file(input_file_name, fformat="roff") + return obs_field_object + + +def read_upscaled_field_from_file(response_settings, iteration): + input_file_name = response_settings["upscaled_file_name"] + if iteration == 0: + filename = "init_files/" + input_file_name + else: + filename = input_file_name + field_object = xtgeo.gridproperty_from_file(filename, fformat="roff") + return field_object + + +def write_obs_pred_diff_field(upscaled_field_object, observation_field_object): + nx, ny, nz = upscaled_field_object.dimensions + values_diff = upscaled_field_object.values - observation_field_object.values + + diff_object = xtgeo.grid3d.GridProperty( + ncol=nx, + nrow=ny, + nlay=nz, + values=values_diff, + discrete=False, + name="DiffObsPred", + ) + + filename = "DiffObsPred.roff" + print( + f"Write field with difference between observation and prediction: {filename} " + ) + diff_object.to_file(filename, fformat="roff") diff --git a/tests/jobs/localisation/example_case/scripts/init_test_case.py b/tests/jobs/localisation/example_case/scripts/init_test_case.py new file mode 100755 index 000000000..07679a090 --- /dev/null +++ b/tests/jobs/localisation/example_case/scripts/init_test_case.py @@ -0,0 +1,283 @@ +#!/usr/bin/env python +""" +Script initialize the test case by creating the grid files, observation files etc +""" +import math +import os +import random + +import xtgeo + + +from common_functions import ( + specify_settings, + write_upscaled_field, + generate_field_and_upscale, +) + + +# pylint: disable=too-many-arguments,invalid-name,missing-function-docstring +# pylint: disable=too-many-locals,import-error + + +def generate_seed_file( + seed_file_name: str = "randomseeds.txt", + start_seed: int = 9828862224, + number_of_seeds: int = 1000, +): + # pylint: disable=unused-variable + print(f"Generate random seed file: {seed_file_name}") + random.seed(start_seed) + with open(seed_file_name, "w", encoding="utf8") as file: + for i in range(number_of_seeds): + file.write(f"{random.randint(1, 999999999)}\n") + + +def obs_positions(grid_size, response_settings, observation_settings): + NX, NY, _ = response_settings["grid_dimension"] + use_eclipse_origin = grid_size["use_eclipse_grid_index_origin"] + + xsize = grid_size["xsize"] + ysize = grid_size["ysize"] + dx = xsize / NX + dy = ysize / NY + cell_indx_list = observation_settings["selected_grid_cells"] + if use_eclipse_origin: + print("Grid index origin: Eclipse standard") + else: + print("Grid index origin: RMS standard") + print( + "Observation reference point coordinates is always " + "from origin at lower left corner" + ) + + pos_list = [] + for indices in cell_indx_list: + Iindx = indices[0] - 1 + Jindx = indices[1] - 1 + x = (Iindx + 0.5) * dx + if use_eclipse_origin: + y = ysize - (Jindx + 0.5) * dy + else: + y = (Jindx + 0.5) * dy + + pos_list.append((x, y)) + + return pos_list + + +def write_localisation_config( + observation_settings, + field_settings, + grid_size, + response_settings, + config_file_name="local_config.yml", + write_scaling=True, +): + obs_index_list = observation_settings["selected_grid_cells"] + field_name = field_settings["field_name"] + corr_ranges = field_settings["field_correlation_range"] + azimuth = field_settings["field_correlation_azimuth"] + space = " " * 2 + space2 = " " * 4 + space3 = " " * 6 + positions = obs_positions(grid_size, response_settings, observation_settings) + print(f"Write localisation config file: {config_file_name}") + with open(config_file_name, "w", encoding="utf8") as file: + file.write("log_level: 3\n") + file.write(f"write_scaling_factors: {write_scaling}\n") + file.write("correlations:\n") + for i, indx in enumerate(obs_index_list): + Iindx, Jindx, Kindx = indx + obs_name = f"OBS_{Iindx}_{Jindx}_{Kindx}" + pos = positions[i] + file.write(f"{space}- name: CORR_{i}\n") + file.write(f"{space2}obs_group:\n") + file.write(f'{space3}add: ["{obs_name}"]\n') + file.write(f"{space2}param_group:\n") + file.write(f'{space3}add: ["{field_name}"]\n') + file.write(f"{space2}field_scale:\n") + file.write(f"{space3}method: gaussian_decay\n") + file.write(f"{space3}main_range: {corr_ranges[0]}\n") + file.write(f"{space3}perp_range: {corr_ranges[1]}\n") + file.write(f"{space3}azimuth: {azimuth}\n") + file.write(f"{space3}ref_point: [ {pos[0]}, {pos[1]} ]\n") + + +def write_gen_obs(upscaled_values, observation_settings): + observation_dir = observation_settings["observation_dir"] + obs_file_name = observation_settings["observation_file"] + obs_data_dir = observation_settings["observation_data_dir"] + cell_indx_list = observation_settings["selected_grid_cells"] + rel_err = observation_settings["rel_error"] + min_err = observation_settings["min_abs_error"] + if not os.path.exists(observation_dir): + print(f"Create directory: {observation_dir} ") + os.makedirs(observation_dir) + data_dir = observation_dir + "/" + obs_data_dir + if not os.path.exists(data_dir): + print(f"Create directory: {data_dir} ") + os.makedirs(data_dir) + + print(f"Write observation file: {obs_file_name} ") + filename = observation_dir + "/" + obs_file_name + with open(filename, "w", encoding="utf8") as obs_file: + number = 0 + for indices in cell_indx_list: + Iindx = indices[0] - 1 + Jindx = indices[1] - 1 + Kindx = indices[2] - 1 + + value = upscaled_values[Iindx, Jindx, Kindx] + value_err = math.fabs(value) * rel_err + value_err = max(value_err, min_err) + + obs_data_relative_file_name = ( + obs_data_dir + + "/obs_" + + str(Iindx + 1) + + "_" + + str(Jindx + 1) + + "_" + + str(Kindx + 1) + + ".txt" + ) + + obs_file.write(f"GENERAL_OBSERVATION OBS_{Iindx+1}_{Jindx+1}_{Kindx+1} ") + obs_file.write("{ ") + obs_file.write( + f"DATA = RESULT_UPSCALED_FIELD ; INDEX_LIST = {number} ; RESTART = 0; " + ) + obs_file.write(f"OBS_FILE = ./{obs_data_relative_file_name} ; ") + obs_file.write(" };\n") + number += 1 + data_file_name = observation_dir + "/" + obs_data_relative_file_name + print(f"Write file: {data_file_name} ") + with open(data_file_name, "w", encoding="utf8") as data_file: + data_file.write(f"{value} {value_err}\n") + + +def create_grid(grid_size, field_settings): + grid_file_name = field_settings["grid_file_name"] + nx, ny, nz = field_settings["grid_dimension"] + xsize = grid_size["xsize"] + ysize = grid_size["ysize"] + zsize = grid_size["zsize"] + if grid_size["use_eclipse_grid_index_origin"]: + flip = -1 + x0 = 0.0 + y0 = ysize + z0 = 0.0 + else: + flip = 1 + x0 = 0.0 + y0 = 0.0 + z0 = 0.0 + + dx = xsize / nx + dy = ysize / ny + dz = zsize / nz + + grid_object = xtgeo.create_box_grid( + dimension=(nx, ny, nz), + origin=(x0, y0, z0), + increment=(dx, dy, dz), + rotation=0.0, + flip=flip, + ) + + print(f"Write grid file: {grid_file_name} ") + grid_object.to_file(grid_file_name, fformat="egrid") + return grid_object + + +def create_upscaled_grid(grid_size, response_settings): + grid_file_name = response_settings["grid_file_name"] + nx, ny, nz = response_settings["grid_dimension"] + xsize = grid_size["xsize"] + ysize = grid_size["ysize"] + zsize = grid_size["zsize"] + if grid_size["use_eclipse_grid_index_origin"]: + flip = -1 + x0 = 0.0 + y0 = ysize + z0 = 0.0 + else: + flip = 1 + x0 = 0.0 + y0 = 0.0 + z0 = 0.0 + + dx = xsize / nx + dy = ysize / ny + dz = zsize / nz + + grid_object = xtgeo.create_box_grid( + dimension=(nx, ny, nz), + origin=(x0, y0, z0), + increment=(dx, dy, dz), + rotation=0.0, + flip=flip, + ) + + print(f"Write grid file: {grid_file_name} ") + grid_object.to_file(grid_file_name, fformat="egrid") + return grid_object + + +def main(): + """ + Initialize seed file, grid files, observation files and localisation config file + """ + # Settings are specified here + ( + grid_size, + field_settings, + response_settings, + observation_settings, + ) = specify_settings() + + # Create seed file + generate_seed_file(field_settings["field_seed_file"]) + + # Create grid for the field parameter + create_grid(grid_size, field_settings) + + # Create coarse grid to be used in QC of upscaled field parameter + create_upscaled_grid(grid_size, response_settings) + + print("Generate field parameter and upscale this.") + print( + f"The upscaled field {observation_settings['3D_param_file_name']} " + "is used when extracting observations." + ) + + # Simulate field (with trend) + real_number = 0 + upscaled_values = generate_field_and_upscale( + field_settings, grid_size, response_settings, observation_settings, real_number + ) + + # Create observations by extracting from existing upscaled field + write_gen_obs(upscaled_values, observation_settings) + + # Write upscaled field used as truth realisation + write_upscaled_field( + upscaled_values, + observation_settings["3D_param_file_name"], + selected_cell_index_list=observation_settings["selected_grid_cells"], + ) + + # Write file for non-adaptive localisation using distance based localisation + write_localisation_config( + observation_settings, + field_settings, + grid_size, + response_settings, + config_file_name="local_config.yml", + write_scaling=True, + ) + + +if __name__ == "__main__": + main() diff --git a/tests/jobs/localisation/example_case/scripts/sim_fields.py b/tests/jobs/localisation/example_case/scripts/sim_fields.py index defcde729..7a44aa305 100755 --- a/tests/jobs/localisation/example_case/scripts/sim_fields.py +++ b/tests/jobs/localisation/example_case/scripts/sim_fields.py @@ -2,15 +2,19 @@ """ Script used as forward model in ERT to test localisation. """ -import math import os -import random import sys -import xtgeo -import gaussianfft as sim -import numpy as np - +# pylint: disable=import-error +from common_functions import ( + specify_settings, + generate_field_and_upscale, + read_field_from_file, + upscaling, + read_obs_field_from_file, + read_upscaled_field_from_file, + write_obs_pred_diff_field, +) # NOTE: xtgeo MUST be imported BEFORE gaussianfft # The reason is that xtgeo has functions importing roxar API and @@ -46,361 +50,7 @@ # it is imported. The most robust solution would be to generate gaussianfft # to use exactly the same. - - -# pylint: disable=missing-function-docstring, too-many-locals, invalid-name -# pylint: disable=bare-except, raise-missing-from -# pylint: disable= redefined-outer-name, too-many-nested-blocks - - -def specify_settings(): - """ - grid_size - length, width, thickness of a box containing the field - Same size is used for both fine scale grid with the simulated field - and the coarse scale grid containing upscaled values - of the simulated field. - - field_settings - Define the dimension (number of grid cells) for fine scale grid, - name of output files and specification of model parameters for - simulation of gaussian field with option to use linear trend. - Relative standard deviation specify standard deviation of - gaussian residual field relative to the trends span of value - (max trend value - min trend value) - - response_settings - Specify the coarse grid dimensions, name of file and type - of average operation to calculated upscaled values that - are predictions of observations of the same grid cells. - Which cell indices are observed are specified in - observation settings. - - observation_settings - Specify name of files for generated observations - and also which grid cells from coarse grid is used - as observables. - (Cells that have values that are used as observations) - """ - grid_size = { - "xsize": 7500.0, - "ysize": 12500.0, - "zsize": 50.0, - "use_eclipse_grid_index_origin": True, - } - field_settings = { - "grid_dimension": [150, 250, 1], - "grid_file_name": "GRID.EGRID", - "field_correlation_range": [3000.0, 1000.0, 2.0], - "field_correlation_azimuth": 45.0, - "field_correlation_dip": 0.0, - "field_variogram": "spherical", - "field_trend_params": [1.0, -1.0], - "field_trend_relstd": 0.05, - "field_trend_use": 0, - "field_name": "FIELDPARAM", - "field_updated_file_name": "FieldParam.roff", - "field_initial_file_name": "init_files/FieldParam.roff", - "field_seed_file": "randomseeds.txt", - } - response_settings = { - "grid_dimension": [15, 25, 1], - "upscaled_file_name": "Upscaled.roff", - "grid_file_name": "UpscaleGrid.EGRID", - "response_function": "average", - "gen_data_file_name": "UpscaledField_0.txt", - "all": True, - } - - observation_settings = { - "observation_dir": "observations", - "observation_file": "observations.obs", - "observation_data_dir": "obs_data", - "3D_param_file_name": "init_files/UpscaledObsField.roff", - "rel_error": 0.10, - "min_abs_error": 0.01, - "selected_grid_cells": [ - [5, 10, 1], - ], - # "selected_grid_cells":[ - # [15, 1, 1], - # [15, 25, 1], - # [ 1, 25, 1], - # [ 1, 1, 1], - # [ 3, 3, 1], - # [ 5, 23, 1], - # [11, 4, 1], - # [ 3, 12, 1], - # [13, 18, 1], - # ], - } - return grid_size, field_settings, response_settings, observation_settings - - -def generate_seed_file( - seed_file_name: str = "randomseeds.txt", - start_seed: int = 9828862224, - number_of_seeds: int = 1000, -): - # pylint: disable=unused-variable - print(f"Generate random seed file: {seed_file_name}") - random.seed(start_seed) - with open(seed_file_name, "w", encoding="utf8") as file: - for i in range(number_of_seeds): - file.write(f"{random.randint(1, 999999999)}\n") - - -def get_seed(seed_file_name, r_number): - with open(seed_file_name, "r", encoding="utf8") as file: - lines = file.readlines() - try: - seed_value = int(lines[r_number - 1]) - except: # noqa: E722 - raise IOError("Seed value not found for realization {r_number} ") - return seed_value - - -def obs_positions(grid_size, response_settings, observation_settings): - NX, NY, _ = response_settings["grid_dimension"] - use_eclipse_origin = grid_size["use_eclipse_grid_index_origin"] - - xsize = grid_size["xsize"] - ysize = grid_size["ysize"] - dx = xsize / NX - dy = ysize / NY - cell_indx_list = observation_settings["selected_grid_cells"] - if use_eclipse_origin: - print("Grid index origin: Eclipse standard") - else: - print("Grid index origin: RMS standard") - print( - "Observation reference point coordinates is always " - "from origin at lower left corner" - ) - - pos_list = [] - for indices in cell_indx_list: - Iindx = indices[0] - 1 - Jindx = indices[1] - 1 - x = (Iindx + 0.5) * dx - if use_eclipse_origin: - y = ysize - (Jindx + 0.5) * dy - else: - y = (Jindx + 0.5) * dy - - pos_list.append((x, y)) - - return pos_list - - -def write_localisation_config( - observation_settings, - field_settings, - positions, - config_file_name="local_config.yml", - write_scaling=True, -): - obs_index_list = observation_settings["selected_grid_cells"] - field_name = field_settings["field_name"] - corr_ranges = field_settings["field_correlation_range"] - azimuth = field_settings["field_correlation_azimuth"] - space = " " * 2 - space2 = " " * 4 - space3 = " " * 6 - print(f"Write localisation config file: {config_file_name}") - with open(config_file_name, "w", encoding="utf8") as file: - file.write("log_level: 3\n") - file.write(f"write_scaling_factors: {write_scaling}\n") - file.write("correlations:\n") - for i, indx in enumerate(obs_index_list): - I, J, K = indx - obs_name = f"OBS_{I}_{J}_{K}" - pos = positions[i] - file.write(f"{space}- name: CORR_{i}\n") - file.write(f"{space2}obs_group:\n") - file.write(f'{space3}add: ["{obs_name}"]\n') - file.write(f"{space2}param_group:\n") - file.write(f'{space3}add: ["{field_name}"]\n') - file.write(f"{space2}field_scale:\n") - file.write(f"{space3}method: gaussian_decay\n") - file.write(f"{space3}main_range: {corr_ranges[0]}\n") - file.write(f"{space3}perp_range: {corr_ranges[1]}\n") - file.write(f"{space3}azimuth: {azimuth}\n") - file.write(f"{space3}ref_point: [ {pos[0]}, {pos[1]} ]\n") - - -def upscaling( - field_values, response_settings, observation_settings, write_field=True, iteration=0 -): - response_function_name = response_settings["response_function"] - upscaled_file_name = response_settings["upscaled_file_name"] - NX, NY, NZ = response_settings["grid_dimension"] - calculate_all = response_settings["all"] - - coarse_cell_index_list = observation_settings["selected_grid_cells"] - upscaled_values = np.zeros((NX, NY, NZ), dtype=np.float32, order="F") - upscaled_values[:, :, :] = -999 - - if response_function_name == "average": - upscaled_values = upscale_average( - field_values, coarse_cell_index_list, upscaled_values, use_all=calculate_all - ) - - if iteration == 0: - upscaled_file_name = "init_files/" + upscaled_file_name - - if write_field: - write_upscaled_field(upscaled_values, upscaled_file_name) - return upscaled_values - - -def write_upscaled_field( - upscaled_values, upscaled_file_name, selected_cell_index_list=None -): - nx, ny, nz = upscaled_values.shape - field_name = "Upscaled" - - field_object = xtgeo.grid3d.GridProperty( - ncol=nx, - nrow=ny, - nlay=nz, - values=upscaled_values, - discrete=False, - name=field_name, - ) - - print(f"Write upscaled field file: {upscaled_file_name} ") - field_object.to_file(upscaled_file_name, fformat="roff") - if selected_cell_index_list is not None: - selected_upscaled_values = np.zeros((nx, ny, nz), dtype=np.float32, order="F") - selected_upscaled_values[:, :, :] = -1 - for indices in selected_cell_index_list: - Iindx = indices[0] - 1 - Jindx = indices[1] - 1 - Kindx = indices[2] - 1 - selected_upscaled_values[Iindx, Jindx, Kindx] = upscaled_values[ - Iindx, Jindx, Kindx - ] - - field_name_selected = field_name + "_conditioned_cells" - file_name_selected = "init_files/" + field_name_selected + ".roff" - cond_field_object = xtgeo.grid3d.GridProperty( - ncol=nx, - nrow=ny, - nlay=nz, - values=selected_upscaled_values, - discrete=False, - name=field_name_selected, - ) - print(f"Write conditioned cell values as field: {file_name_selected}") - cond_field_object.to_file(file_name_selected, fformat="roff") - - return field_object - - -def upscale_average( - field_values, coarse_cell_index_list, upscaled_values, use_all=False -): - """ - Input: field_values (numpy 3D) - coarse_cell_index_list (list of tuples (I,J,K)) - Output: upscaled_values (numpy 3D) initialized outside - but filled in specified (I,J,K) cells. - """ - nx, ny, nz = field_values.shape - NX, NY, NZ = upscaled_values.shape - - print(f"Number of fine scale grid cells: (nx,ny,nz): ({nx},{ny},{nz})") - print(f"Number of coarse scale grid cells: (NX,NY,NZ): ({NX},{NY},{NZ}) ") - mx = int(nx / NX) - my = int(ny / NY) - mz = int(nz / NZ) - print( - "Number of fine scale grid cells per coarse grid cell: " - f"(mx,my,mz): ({mx},{my},{mz}) " - ) - if use_all: - print("Calculate upscaled values for all grid cells") - for Kindx in range(NZ): - for Jindx in range(NY): - for Iindx in range(NX): - istart = mx * Iindx - iend = istart + mx - jstart = my * Jindx - jend = jstart + my - kstart = mz * Kindx - kend = kstart + mz - sum_val = 0.0 - for k in range(kstart, kend): - for j in range(jstart, jend): - for i in range(istart, iend): - sum_val += field_values[i, j, k] - upscaled_values[Iindx, Jindx, Kindx] = sum_val / (mx * my * mz) - - else: - print("Calculate upscaled values for selected grid cells") - for indices in coarse_cell_index_list: - Iindx = indices[0] - 1 - Jindx = indices[1] - 1 - Kindx = indices[2] - 1 - istart = mx * Iindx - iend = istart + mx - jstart = my * Jindx - jend = jstart + my - kstart = mz * Kindx - kend = kstart + mz - sum_val = 0.0 - for k in range(kstart, kend): - for j in range(jstart, jend): - for i in range(istart, iend): - sum_val += field_values[i, j, k] - upscaled_values[Iindx, Jindx, Kindx] = sum_val / (mx * my * mz) - - return upscaled_values - - -def write_gen_obs(upscaled_values, observation_settings): - observation_dir = observation_settings["observation_dir"] - obs_file_name = observation_settings["observation_file"] - obs_data_dir = observation_settings["observation_data_dir"] - cell_indx_list = observation_settings["selected_grid_cells"] - rel_err = observation_settings["rel_error"] - min_err = observation_settings["min_abs_error"] - print(f"Write observation file: {obs_file_name} ") - if not os.path.exists(obs_data_dir): - print(f"Create directory: {obs_data_dir} ") - os.makedirs(obs_data_dir) - filename = observation_dir + "/" + obs_file_name - with open(filename, "w", encoding="utf8") as obs_file: - number = 0 - for indices in cell_indx_list: - Iindx = indices[0] - 1 - Jindx = indices[1] - 1 - Kindx = indices[2] - 1 - - value = upscaled_values[Iindx, Jindx, Kindx] - value_err = math.fabs(value) * rel_err - value_err = max(value_err, min_err) - - obs_data_relative_file_name = ( - obs_data_dir - + "/obs_" - + str(Iindx + 1) - + "_" - + str(Jindx + 1) - + "_" - + str(Kindx + 1) - + ".txt" - ) - - obs_file.write(f"GENERAL_OBSERVATION OBS_{Iindx+1}_{Jindx+1}_{Kindx+1} ") - obs_file.write("{ ") - obs_file.write( - f"DATA = RESULT_UPSCALED_FIELD ; INDEX_LIST = {number} ; RESTART = 0; " - ) - obs_file.write(f"OBS_FILE = ./{obs_data_relative_file_name} ; ") - obs_file.write(" };\n") - number += 1 - data_file_name = observation_dir + "/" + obs_data_relative_file_name - print(f"Write file: {data_file_name} ") - with open(data_file_name, "w", encoding="utf8") as data_file: - data_file.write(f"{value} {value_err}\n") +# pylint: disable=missing-function-docstring,invalid-name def write_prediction_gen_data(upscaled_values, observation_settings, response_settings): @@ -418,212 +68,40 @@ def write_prediction_gen_data(upscaled_values, observation_settings, response_se file.write(f"{value}\n") -def trend(grid_size, field_settings): - """ - Return 3D numpy array with values following a linear trend - scaled to take values between 0 and 1. - """ - nx, ny, nz = field_settings["grid_dimension"] - xsize = grid_size["xsize"] - ysize = grid_size["ysize"] - a, b = field_settings["field_trend_params"] - - x0 = 0.0 - y0 = 0.0 - dx = xsize / nx - dy = ysize / ny - - maxsize = ysize - if xsize > ysize: - maxsize = xsize - - val = np.zeros((nx, ny, nz), dtype=np.float32, order="F") - for i in range(nx): - x = x0 + i * dx - for j in range(ny): - y = y0 + j * dy - for k in range(nz): - val[i, j, k] = a * (x - x0) / maxsize + b * (y - y0) / maxsize - - minval = np.min(val) - maxval = np.max(val) - val_normalized = (val - minval) / (maxval - minval) - return val_normalized - - -def simulate_field(grid_size, field_settings, start_seed): - # pylint: disable=no-member, - variogram_name = field_settings["field_variogram"] - corr_ranges = field_settings["field_correlation_range"] - xrange = corr_ranges[0] - yrange = corr_ranges[1] - zrange = corr_ranges[2] - - azimuth = field_settings["field_correlation_azimuth"] - dip = field_settings["field_correlation_dip"] - - nx, ny, nz = field_settings["grid_dimension"] - xsize = grid_size["xsize"] - ysize = grid_size["ysize"] - zsize = grid_size["zsize"] - - dx = xsize / nx - dy = ysize / ny - dz = zsize / nz - - print(f"Start seed: {start_seed}") - sim.seed(start_seed) - - variogram = sim.variogram( - variogram_name, - xrange, - perp_range=yrange, - depth_range=zrange, - azimuth=azimuth - 90, - dip=dip, - ) - - print(f"Simulate field with size: nx={nx},ny={ny} ") - field1D = sim.simulate(variogram, nx, dx, ny, dy, nz, dz) - field = field1D.reshape((nx, ny, nz), order="F") - return field - - -def create_grid(grid_size, field_settings): - grid_file_name = field_settings["grid_file_name"] - nx, ny, nz = field_settings["grid_dimension"] - xsize = grid_size["xsize"] - ysize = grid_size["ysize"] - zsize = grid_size["zsize"] - if grid_size["use_eclipse_grid_index_origin"]: - flip = -1 - x0 = 0.0 - y0 = ysize - z0 = 0.0 +def get_iteration_and_real_number(argv): + iter_env_name = "_ERT_ITERATION_NUMBER" + if iter_env_name in os.environ: + iteration = int(os.environ.get(iter_env_name)) else: - flip = 1 - x0 = 0.0 - y0 = 0.0 - z0 = 0.0 - - dx = xsize / nx - dy = ysize / ny - dz = zsize / nz - - grid_object = xtgeo.create_box_grid( - dimension=(nx, ny, nz), - origin=(x0, y0, z0), - increment=(dx, dy, dz), - rotation=0.0, - flip=flip, - ) - - print(f"Write grid file: {grid_file_name} ") - grid_object.to_file(grid_file_name, fformat="egrid") - return grid_object - - -def create_upscaled_grid(grid_size, response_settings): - grid_file_name = response_settings["grid_file_name"] - nx, ny, nz = response_settings["grid_dimension"] - xsize = grid_size["xsize"] - ysize = grid_size["ysize"] - zsize = grid_size["zsize"] - if grid_size["use_eclipse_grid_index_origin"]: - flip = -1 - x0 = 0.0 - y0 = ysize - z0 = 0.0 - else: - flip = 1 - x0 = 0.0 - y0 = 0.0 - z0 = 0.0 - - dx = xsize / nx - dy = ysize / ny - dz = zsize / nz - - grid_object = xtgeo.create_box_grid( - dimension=(nx, ny, nz), - origin=(x0, y0, z0), - increment=(dx, dy, dz), - rotation=0.0, - flip=flip, - ) - - print(f"Write grid file: {grid_file_name} ") - grid_object.to_file(grid_file_name, fformat="egrid") - return grid_object - - -def export_field(field_settings, field3D): - # Export initial ensemble field - nx, ny, nz = field_settings["grid_dimension"] - field_name = field_settings["field_name"] - field_file_name = field_settings["field_initial_file_name"] - - field_object = xtgeo.grid3d.GridProperty( - ncol=nx, nrow=ny, nlay=nz, values=field3D, discrete=False, name=field_name - ) - - print(f"Write field file: {field_file_name} ") - field_object.to_file(field_file_name, fformat="roff") - return field_object - - -def read_field_from_file(field_settings): - input_file_name = field_settings["field_updated_file_name"] - name = field_settings["field_name"] - field_object = xtgeo.gridproperty_from_file( - input_file_name, fformat="roff", name=name - ) - return field_object - - -def read_obs_field_from_file(observation_settings): - input_file_name = observation_settings["3D_param_file_name"] - obs_field_object = xtgeo.gridproperty_from_file(input_file_name, fformat="roff") - return obs_field_object - + if len(argv) < 3: + raise IOError("Missing command line arguments ") + arg1 = argv[1] + if arg1 is None: + raise IOError( + "Missing iteration number (argv[1]) when running this script manually" + ) + iteration = int(arg1) + print(f"ERT iteration: {iteration}") -def read_upscaled_field_from_file(response_settings, iteration): - input_file_name = response_settings["upscaled_file_name"] - if iteration == 0: - filename = "init_files/" + input_file_name + real_env_name = "_ERT_REALIZATION_NUMBER" + if real_env_name in os.environ: + real_number = int(os.environ.get(real_env_name)) else: - filename = input_file_name - field_object = xtgeo.gridproperty_from_file(filename, fformat="roff") - return field_object - - -def write_obs_pred_diff_field(upscaled_field_object, observation_field_object): - nx, ny, nz = upscaled_field_object.dimensions - values_diff = upscaled_field_object.values - observation_field_object.values - - diff_object = xtgeo.grid3d.GridProperty( - ncol=nx, - nrow=ny, - nlay=nz, - values=values_diff, - discrete=False, - name="DiffObsPred", - ) - - filename = "DiffObsPred.roff" - print( - f"Write field with difference between observation and prediction: {filename} " - ) - diff_object.to_file(filename, fformat="roff") + if len(argv) < 3: + raise IOError("Missing command line arguments ") + arg2 = argv[2] + if arg2 is None: + raise IOError( + "Missing real_number (argv[2]) when running this script manually" + ) + real_number = int(arg2) + print(f"ERT realization: {real_number}") + return iteration, real_number def main( - iteration, - real_number, - write_fine_grid=False, - write_coarse_grid=False, + argv, write_upscaled_to_file=False, - extract_and_write_obs=False, write_obs_pred_diff_field_file=False, ): # pylint: disable=too-many-arguments @@ -654,39 +132,16 @@ def main( observation_settings, ) = specify_settings() - # Create and write grid file for fine scale grid - if write_fine_grid: - create_grid(grid_size, field_settings) - - if write_coarse_grid: - create_upscaled_grid(grid_size, response_settings) - + iteration, real_number = get_iteration_and_real_number(argv) if iteration == 0: print(f"Generate new field parameter realization:{real_number} ") # Simulate field (with trend) - seed_file_name = field_settings["field_seed_file"] - relative_std = field_settings["field_trend_relstd"] - start_seed = get_seed(seed_file_name, real_number) - residual_field = simulate_field(grid_size, field_settings, start_seed) - trend_field = trend(grid_size, field_settings) - use_trend = field_settings["field_trend_use"] - - if use_trend == 1: - field3D = trend_field + relative_std * residual_field - else: - field3D = residual_field - - # Write field parameter for fine scale grid - field_object = export_field(field_settings, field3D) - field_values = field_object.values - - # Calculate upscaled values for selected coarse grid cells - upscaled_values = upscaling( - field_values, + upscaled_values = generate_field_and_upscale( + field_settings, + grid_size, response_settings, observation_settings, - write_field=write_upscaled_to_file, - iteration=iteration, + real_number, ) else: @@ -702,27 +157,10 @@ def main( write_field=write_upscaled_to_file, iteration=iteration, ) + # Write GEN_DATA file + write_prediction_gen_data(upscaled_values, observation_settings, response_settings) - if extract_and_write_obs: - # ERT obs files - write_gen_obs(upscaled_values, observation_settings) - - # Write upscaled field used as truth realisation - write_upscaled_field( - upscaled_values, - observation_settings["3D_param_file_name"], - selected_cell_index_list=observation_settings["selected_grid_cells"], - ) - - # Write file for non-adaptive localisation using distance based localisation - positions = obs_positions(grid_size, response_settings, observation_settings) - write_localisation_config( - observation_settings, - field_settings, - positions, - config_file_name="local_config.yml", - ) - + # Optional output if write_obs_pred_diff_field_file: obs_field_object = read_obs_field_from_file(observation_settings) upscaled_field_object = read_upscaled_field_from_file( @@ -730,50 +168,15 @@ def main( ) write_obs_pred_diff_field(upscaled_field_object, obs_field_object) - write_prediction_gen_data(upscaled_values, observation_settings, response_settings) - if __name__ == "__main__": - # Create file with one seed pr realization - make_seed_file = False - - # Create and write fine scale grid file to be used in ERT GRID keyword - write_fine_grid = False - - # Create upscaled grid file to be used in visualization of upscaled values in RMS - write_coarse_grid = False - - # Write 3D parameter of upscaled values for the selected + # Optional output: Write 3D parameter of upscaled values for the selected # grid cells having observations - write_upscaled_to_file = True - write_obs_pred_diff_field_file = True - - # Extract and write observation files for ERT. - # Observations extracted for a 3D parameter from the coarse grid. - extract_and_write_obs = False - - try: - iteration = int(os.environ.get("_ERT_ITERATION_NUMBER")) - extract_and_write_obs = False - except: # noqa: E722 - iteration = int(sys.argv[1]) - print(f"ERT iteration: {iteration}") - - try: - real_number = int(os.environ.get("_ERT_REALIZATION_NUMBER")) - except: # noqa: E722 - real_number = int(sys.argv[2]) - print(f"ERT realization: {real_number}") - - if make_seed_file: - generate_seed_file() + define_write_upscaled_to_file = True + define_write_obs_pred_diff_field_file = True main( - iteration, - real_number, - write_fine_grid=write_fine_grid, - write_coarse_grid=write_coarse_grid, - write_upscaled_to_file=write_upscaled_to_file, - extract_and_write_obs=extract_and_write_obs, - write_obs_pred_diff_field_file=write_obs_pred_diff_field_file, + sys.argv, + write_upscaled_to_file=define_write_upscaled_to_file, + write_obs_pred_diff_field_file=define_write_obs_pred_diff_field_file, )