Modified test case for non-adaptive localisation

semeio/workflows/localisation/local_script_lib.py

@@ -17,6 +17,8 @@
 from ecl.util.geometry import Surface
 from ert.analysis.row_scaling import RowScaling
 from ert.config import Field, GenDataConfig, GenKwConfig, SurfaceConfig
+from ert.field_utils import save_field
+from ert.field_utils.field_file_format import FieldFileFormat
 from numpy import ma
 
 from semeio.workflows.localisation.localisation_debug_settings import (
@@ -609,7 +611,9 @@ def add_ministeps(
                             tapering_range = (
                                 corr_spec.field_scale.normalised_tapering_range
                             )
-                        check_if_ref_point_in_grid(ref_pos, grid_for_field)
+                        check_if_ref_point_in_grid(
+                            ref_pos, grid_for_field, log_level=user_config.log_level
+                        )
                         param_for_field = apply_decay(
                             corr_spec.field_scale.method,
                             row_scaling,
@@ -787,15 +791,21 @@ def add_ministeps(
     return update_steps
 
 
-def check_if_ref_point_in_grid(ref_point, grid):
+def check_if_ref_point_in_grid(ref_point, grid, log_level):
     try:
-        grid.find_cell_xy(ref_point[0], ref_point[1], 0)
+        (i_indx, j_indx) = grid.find_cell_xy(ref_point[0], ref_point[1], 0)
     except ValueError as err:
         raise ValueError(
             f"Reference point {ref_point} corresponds to undefined grid cell "
             f"or is outside the area defined by the grid {grid.get_name()}\n"
-            "Check specification of reference point."
+            "Check specification of reference point ",
+            "and grid index origin of grid with field parameters.",
         ) from err
+    debug_print(
+        f"Reference point {ref_point} has grid indices: ({i_indx}, {j_indx})",
+        LogLevel.LEVEL3,
+        log_level,
+    )
 
 
 @dataclass
@@ -924,35 +934,56 @@ def write_qc_parameter(
         grid,
         param_for_field,
         log_level=LogLevel.OFF,
+        file_format=FieldFileFormat.GRDECL,
     ):
         # pylint: disable=too-many-arguments
+
         if param_for_field is None or field_scale is None:
             return
 
-        scaling_values = np.reshape(
-            param_for_field, (grid.getNX(), grid.getNY(), grid.getNZ()), "F"
+        # Write scaling parameter  once per corr_name
+        if corr_name == cls.corr_name:
+            return
+
+        cls.corr_name = corr_name
+
+        # Need a parameter name <= 8 character long for GRDECL format
+        scaling_kw_name = "S_" + str(cls.scaling_param_number)
+        file_name_without_suffix = (
+            cls.corr_name + "_" + node_name + "_" + scaling_kw_name
         )
 
-        # Write scaling parameter  once per corr_name
-        if corr_name != cls.corr_name:
-            cls.corr_name = corr_name
-            # Need a parameter name <= 8 character long
-            scaling_kw_name = "S_" + str(cls.scaling_param_number)
-            scaling_kw = grid.create_kw(scaling_values, scaling_kw_name, False)
-            filename = (
-                cls.corr_name + "_" + node_name + "_" + scaling_kw_name + ".GRDECL"
-            )
-            print(
-                "Write calculated scaling factor  with name: "
-                f"{scaling_kw_name} to file: {filename}"
-            )
-            debug_print(
-                f"Write calculated scaling factor with name: "
-                f"{scaling_kw_name} to file: {filename}",
-                LogLevel.LEVEL3,
-                log_level,
+        if file_format.upper() == FieldFileFormat.GRDECL:
+            scaling_values = np.reshape(
+                param_for_field, (grid.getNX(), grid.getNY(), grid.getNZ()), "F"
             )
+
+            scaling_kw = grid.create_kw(scaling_values, scaling_kw_name, False)
+            filename = file_name_without_suffix + ".GRDECL"
+
             with cwrap.open(filename, "w") as file:
                 grid.write_grdecl(scaling_kw, file)
-            # Increase parameter number to define unique parameter name
-            cls.scaling_param_number = cls.scaling_param_number + 1
+            # For testing
+            name = scaling_kw_name + "_use_save_field"
+            save_field(scaling_values, name, filename, FieldFileFormat.GRDECL)
+
+        elif file_format.upper() == FieldFileFormat.ROFF:
+            scaling_values = np.reshape(
+                param_for_field, (grid.getNX(), grid.getNY(), grid.getNZ()), "C"
+            )
+            filename = file_name_without_suffix + ".roff"
+            save_field(scaling_values, scaling_kw_name, filename, FieldFileFormat.ROFF)
+
+        print(
+            "Write calculated scaling factor  with name: "
+            f"{scaling_kw_name} to file: {filename}"
+        )
+        debug_print(
+            f"Write calculated scaling factor with name: "
+            f"{scaling_kw_name} to file: {filename}",
+            LogLevel.LEVEL3,
+            log_level,
+        )
+
+        # Increase parameter number to define unique parameter name
+        cls.scaling_param_number = cls.scaling_param_number + 1

tests/jobs/localisation/example_case/README

@@ -27,9 +27,8 @@
 ​
 1. **Preparation:** Prepare ERT config input by running `scripts/init_test_case.py`.
 2. Directories for observations are created automatically according to the default settings.
-3. Make the directory `init_files` if not existing.
-4. Activate or deactivate localisation in ERT config file (`HOOK_WORKFLOW LOAD_WORKFLOW` for localisation).
-5. Now ready to run ERT.
+3. Activate or deactivate localisation in ERT config file (`HOOK_WORKFLOW LOAD_WORKFLOW` for localisation).
+4. Now ready to run ERT.
 ​
 ### What the Script `sim_fields.py` Does:
 ​

tests/jobs/localisation/example_case/scripts/common_functions.py

@@ -26,6 +26,8 @@ class GridSize:
 
     xsize: float = 7500.0
     ysize: float = 12500.0
+    #    xsize: float = 7500.0
+    #    ysize: float = 7500.0
     zsize: float = 50.0
     use_eclipse_grid_index_origo: bool = True
 
@@ -42,13 +44,14 @@ class Field:
     """
 
     # pylint: disable=too-many-instance-attributes
-    name: str = "FIELDPARAM"
+    name: str = "FIELDPAR"
     algorithm: str = "gstools"
-    initial_file_name: str = "init_files/FieldParam.roff"
-    updated_file_name: str = "FieldParam.roff"
+    file_format: str = "GRDECL"
+    initial_file_name: str = "init_files/FieldParam.grdecl"
+    updated_file_name: str = "FieldParam.grdecl"
     seed_file: str = "randomseeds.txt"
     variogram: str = "exponential"
-    correlation_range: Tuple[float] = (5000.0, 3000.0, 2.0)
+    correlation_range: Tuple[float] = (3000.0, 1000.0, 2.0)
     correlation_azimuth: float = 45.0
     correlation_dip: float = 0.0
     correlation_exponent: float = 1.9
@@ -91,7 +94,18 @@ class Observation:
     param_file_name: str = "init_files/UpscaledObsField.roff"
     rel_error: float = 0.10
     min_abs_error: float = 0.01
-    selected_grid_cells: Tuple[Tuple[int]] = ((5, 10, 1), (10, 20, 1))
+    # selected_grid_cells: Tuple[Tuple[int]] = ((5, 10, 1), (10, 5, 1))
+    selected_grid_cells: Tuple[Tuple[int]] = (8, 12, 1)
+
+
+@dataclass
+class Localisation:
+    """
+    Specify settings for the localisation config file.
+    """
+
+    method: str = "gaussian"
+    # method: str = "constant"
 
 
 @dataclass
@@ -115,13 +129,14 @@ class Settings:
     field: Field = Field()
     response: Response = Response()
     observation: Observation = Observation()
+    localisation: Localisation = Localisation()
     optional: Optional = Optional()
 
 
 settings = Settings()
 
 
-def generate_field_and_upscale(real_number):
+def generate_field_and_upscale(real_number, grid):
     seed_file_name = settings.field.seed_file
     relative_std = settings.field.trend_relstd
     use_trend = settings.field.trend_use
@@ -141,7 +156,7 @@ def generate_field_and_upscale(real_number):
         field3D = residual_field
 
     # Write field parameter for fine scale grid
-    field_object = export_field(field3D)
+    field_object = export_field(field3D, grid)
     field_values = field_object.values
 
     # Calculate upscaled values for selected coarse grid cells
@@ -210,14 +225,16 @@ def write_upscaled_field(
     )
 
     print(f"Write upscaled field file: {upscaled_file_name}  ")
+    # Grid index order to xtgeo must be c-order
     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 = np.zeros((nx, ny, nz), dtype=np.float32)
         selected_upscaled_values[:, :, :] = -1
-        for indices in selected_cell_index_list:
-            Iindx = indices[0] - 1
-            Jindx = indices[1] - 1
-            Kindx = indices[2] - 1
+        nobs = get_nobs()
+        for obs_number in range(nobs):
+            (Iindx, Jindx, Kindx) = get_cell_indices(
+                obs_number, nobs, selected_cell_index_list
+            )
             selected_upscaled_values[Iindx, Jindx, Kindx] = upscaled_values[
                 Iindx, Jindx, Kindx
             ]
@@ -371,9 +388,19 @@ def simulate_field(start_seed):
     )
 
     print(f"Simulate field with size: nx={nx},ny={ny} ")
+    # gaussianfft.simulate  will save the values in F-order
     field1D = sim.simulate(variogram, nx, dx, ny, dy, nz, dz)
-    field = field1D.reshape((nx, ny, nz), order="F")
-    return field
+    field_sim = field1D.reshape((nx, ny, nz), order="F")
+    if settings.grid_size.use_eclipse_grid_index_origo:
+        field_c_order = np.zeros((nx, ny, nz), dtype=np.float32)
+        j_indices = -np.arange(ny) + ny - 1
+        # Flip j index and use c-order
+        field_c_order[:, j_indices, :] = field_sim[:, :, :]
+        return field_c_order
+    # Change to C-order
+    field_c_order = np.zeros((nx, ny, nz), dtype=np.float32)
+    field_c_order[:, :, :] = field_sim[:, :, :]
+    return field_c_order
 
 
 def simulate_field_using_gstools(start_seed):
@@ -434,44 +461,71 @@ def simulate_field_using_gstools(start_seed):
     # print(f"Field type name: {srf.name}  ")
     # print(f"Field nugget: {srf.nugget} ")
     # print(f"Field opt arg:   {srf.opt_arg}")
-    field = field_srf.reshape((nx, ny, nz), order="F")
+    field = field_srf.reshape((nx, ny, nz), dtype=np.float32)
     if settings.grid_size.use_eclipse_grid_index_origo:
-        field_result = np.zeros((nx, ny, nz), dtype=np.float32)
+        field_flip_j_index = np.zeros((nx, ny, nz), dtype=np.float32)
         j_indices = -np.arange(ny) + ny - 1
-        field_result[:, j_indices, :] = field[:, :, :]
-        return field_result
+        field_flip_j_index[:, j_indices, :] = field[:, :, :]
+        return field_flip_j_index
 
     return field
 
 
-def export_field(field3D):
+def read_grid():
+    file_format = settings.field.file_format
+    filename = settings.field.grid_file_name
+    if file_format.upper() == "GRDECL":
+        grid = xtgeo.grid_from_file(filename, fformat="egrid")
+    elif file_format.upper() == "ROFF":
+        grid = xtgeo.grid_from_file(filename, fformat="roff")
+    else:
+        raise IOError(f"Unknown grid file format:  {file_format}")
+    return grid
+
+
+def export_field(field3D, grid):
     """
     Export initial realization of field to roff format
+    Input field3D should be C-index order since xtgeo requires that
     """
 
-    nx, ny, nz = settings.field.grid_dimension
+    # nx, ny, nz = settings.field.grid_dimension
     field_name = settings.field.name
     field_file_name = settings.field.initial_file_name
+    file_format = settings.field.file_format.upper()
 
     field_object = xtgeo.grid3d.GridProperty(
-        ncol=nx, nrow=ny, nlay=nz, values=field3D, discrete=False, name=field_name
+        grid, values=field3D, discrete=False, name=field_name
     )
 
     print(f"Write field file: {field_file_name}  ")
-    field_object.to_file(field_file_name, fformat="roff")
+    if file_format == "GRDECL":
+        field_object.to_file(field_file_name, fformat="grdecl", dtype="float32")
+    elif file_format == "ROFF":
+        field_object.to_file(field_file_name, fformat="roff")
+    else:
+        raise IOError(f"Unknown file format for fields: {file_format} ")
     return field_object
 
 
-def read_field_from_file():
+def read_field_from_file(grid):
     """
     Read field from roff formatted file and return xtgeo property object
     """
 
     input_file_name = settings.field.updated_file_name
     name = settings.field.name
-    field_object = xtgeo.gridproperty_from_file(
-        input_file_name, fformat="roff", name=name
-    )
+    file_format = settings.field.file_format.upper()
+    if file_format == "GRDECL":
+        field_object = xtgeo.grid3d.GridProperty(
+            input_file_name, fformat="grdecl", grid=grid, name=name
+        )
+    elif file_format == "ROFF":
+        field_object = xtgeo.gridproperty_from_file(
+            input_file_name, fformat="roff", name=name
+        )
+    else:
+        raise IOError(f"Unknown file format for fields: {file_format} ")
     return field_object
 
 
@@ -524,3 +578,33 @@ def write_obs_pred_diff_field(upscaled_field_object, observation_field_object):
         f"Write field with difference between observation and prediction: {filename}  "
     )
     diff_object.to_file(filename, fformat="roff")
+
+
+def get_cell_indices(obs_number, nobs, cell_indx_list):
+    if nobs == 1:
+        Iindx = cell_indx_list[0] - 1
+        Jindx = cell_indx_list[1] - 1
+        Kindx = cell_indx_list[2] - 1
+    else:
+        Iindx = cell_indx_list[obs_number][0] - 1
+        Jindx = cell_indx_list[obs_number][1] - 1
+        Kindx = cell_indx_list[obs_number][2] - 1
+
+    return (Iindx, Jindx, Kindx)
+
+
+def get_nobs():
+    """
+    Check if cell_indx_list is a single tuple (i,j,k) or
+    a list of tuples  of type (i,j,k).
+    Return number of cell indices found in list
+    """
+    cell_indx_list = settings.observation.selected_grid_cells
+    is_list_of_ints = all(isinstance(indx, int) for indx in cell_indx_list)
+    if is_list_of_ints:
+        nobs = 1
+    else:
+        # list of tuples (i,j,k)
+        nobs = len(cell_indx_list)
+        print(f"nobs: {nobs}  ")
+    return nobs

tests/jobs/localisation/example_case/scripts/import_field_parameters.py

@@ -14,5 +14,9 @@
 FIELD_NAMES = [
     "FieldParam",
 ]
+FILE_FORMAT = "ROFF"
+ITERATION = 3
 
-import_from_scratch_directory(PRJ, GRID_MODEL_NAME, FIELD_NAMES, CASE_NAME, SCRATCH)
+import_from_scratch_directory(
+    PRJ, GRID_MODEL_NAME, FIELD_NAMES, CASE_NAME, SCRATCH, FILE_FORMAT, ITERATION
+)