update

tests/jobs/localisation/example_case/README

@@ -35,11 +35,7 @@ Typical workflow:
     Edit the input yml file to change settings.
 2. Directories for observations are created automatically according to the default (or modified) settings.
 3. Make the directory init_files if not existing.
-4. Run the script init_test_case.py   located under scripts from the example_case dirctory.
-    Run it with or without one argument which is a specification of settings.
-    Default settings are used for settings not specified. The file can be ommited if one want to use the default
-    settings for all parameters. Example file containing all possible settings is the file example_case.yml
-5. If non-default settings are specified and the init_test_case.py is run with an yml file as argument, also the ERT forward
+4. If non-default settings are specified and the init_test_case.py is run with an yml file as argument, also the ERT forward
     model SIM_FIELD need the same input as the third argument for the script sim_fields.py
     The first two arguments are iteration and realisation number. Edit FM_SIM_FIELD and specify wanted settings file (yml file).
 5. Activate/not activate localisation in ERT config file (HOOK_WORKFLOW  LOAD_WORKFLOW for localisation)

tests/jobs/localisation/example_case/example_case.yml

@@ -7,13 +7,15 @@ settings:
 
   field:
     name: "FIELDPARAM"
+    algorithm: "gstools"
     initial_file_name: "init_files/FieldParam.roff"
     updated_file_name: "FieldParam.roff"
     seed_file: "randomseeds.txt"
-    variogram: "gaussian"
-    correlation_range: [5000.0, 3000.0,  2.0]
+    variogram: "exponential"
+    correlation_range: [1000.0, 1000.0,  2.0]
     correlation_azimuth:  45.0
     correlation_dip: 0.0
+    correlation_exponent: 2.0
     trend_use: False
     trend_params: [ 1.0, -1.0 ]
     trend_relstd: 0.15

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

@@ -2,48 +2,13 @@
 Common functions used by the scripts: init_test_case.py and sim_field.py
 """
 import copy
+import math
 
+import gstools as gs
 import numpy as np
 import yaml
 
 import xtgeo  # isort: skip
-import gaussianfft as sim  # isort: skip
-
-
-# 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
@@ -87,13 +52,15 @@ def specify_settings(spec_dict=None, yml_file_name=None):
         },
         "field": {
             "name": "FIELDPARAM",
+            "algorithm": "gstools",
             "initial_file_name": "init_files/FieldParam.roff",
             "updated_file_name": "FieldParam.roff",
             "seed_file": "randomseeds.txt",
-            "variogram": "gaussian",
+            "variogram": "exponential",
             "correlation_range": [5000.0, 3000.0, 2.0],
             "correlation_azimuth": 45.0,
             "correlation_dip": 0.0,
+            "correlation_exponent": 1.9,
             "trend_use": 0,
             "trend_params": [1.0, -1.0],
             "trend_relstd": 0.05,
@@ -167,19 +134,6 @@ def specify_settings(spec_dict=None, yml_file_name=None):
     return settings
 
 
-# def update_settings(settings_dict, spec_dict, main_key, sub_keys):
-#     if main_key in spec_dict:
-#         main_key_dict = spec_dict[main_key]
-#         print(f"main_key_dict:  {main_key_dict}")
-#         if main_key_dict is not None:
-#             print(f"sub_keys:  {sub_keys}")
-#             for key in sub_keys:
-#                 if key in main_key_dict:
-#                     print(f"Settings updated for:{main_key} with sub key:  {key}  ")
-#                     settings_dict[key] = main_key_dict[key]
-#     return settings_dict
-
-
 def read_test_config(config_file_name):
     if config_file_name is None:
         spec_dict = None
@@ -200,12 +154,18 @@ def generate_field_and_upscale(settings, real_number):
     seed_file_name = settings["field"]["seed_file"]
     relative_std = settings["field"]["trend_relstd"]
     use_trend = settings["field"]["trend_use"]
-
+    algorithm_method = settings["field"]["algorithm"]
     start_seed = get_seed(seed_file_name, real_number)
-    residual_field = simulate_field(settings, start_seed)
+    if algorithm_method == "gstools":
+        print(f"Use algorithm: {algorithm_method}")
+        residual_field = simulate_field_using_gstools(settings, start_seed)
+    else:
+        print("Use algorithm: gaussianfft")
+        residual_field = simulate_field(settings, start_seed)
     if use_trend == 1:
         trend_field = trend(settings)
         field3D = trend_field + relative_std * residual_field
+
     else:
         field3D = residual_field
 
@@ -410,7 +370,9 @@ def trend(settings):
 
 
 def simulate_field(settings, start_seed):
-    # pylint: disable=no-member,
+    # pylint: disable=no-member,import-outside-toplevel
+    import gaussianfft as sim  # isort: skip
+
     variogram_name = settings["field"]["variogram"]
     corr_ranges = settings["field"]["correlation_range"]
     xrange = corr_ranges[0]
@@ -419,7 +381,7 @@ def simulate_field(settings, start_seed):
 
     azimuth = settings["field"]["correlation_azimuth"]
     dip = settings["field"]["correlation_dip"]
-
+    alpha = settings["field"]["correlation_exponent"]
     nx, ny, nz = settings["field"]["grid_dimension"]
     xsize = settings["grid_size"]["xsize"]
     ysize = settings["grid_size"]["ysize"]
@@ -439,6 +401,7 @@ def simulate_field(settings, start_seed):
         depth_range=zrange,
         azimuth=azimuth - 90,
         dip=dip,
+        power=alpha,
     )
 
     print(f"Simulate field with size: nx={nx},ny={ny} ")
@@ -447,6 +410,84 @@ def simulate_field(settings, start_seed):
     return field
 
 
+def simulate_field_using_gstools(settings, start_seed):
+    # pylint: disable=no-member,
+
+    variogram_name = settings["field"]["variogram"]
+    corr_ranges = settings["field"]["correlation_range"]
+    xrange = corr_ranges[0]
+    yrange = corr_ranges[1]
+    zrange = corr_ranges[2]
+
+    azimuth = settings["field"]["correlation_azimuth"]
+    alpha = settings["field"]["correlation_exponent"]
+
+    nx, ny, nz = settings["field"]["grid_dimension"]
+    xsize = settings["grid_size"]["xsize"]
+    ysize = settings["grid_size"]["ysize"]
+    zsize = settings["grid_size"]["zsize"]
+
+    dx = xsize / nx
+    dy = ysize / ny
+    dz = zsize / nz
+
+    x = np.arange(0.5 * dx, xsize, dx)
+    y = np.arange(0.5 * dy, ysize, dy)
+    z = np.arange(0.5 * dz, zsize, dz)
+    # Rescale factor is set to:
+    #    sqrt(3.0) for gaussian correlation functions,
+    #    3.0 for exponetial correlation function,
+    #    pow(3.0, 1/alpha) for general exponential correlation function
+    #    with exponent alpha
+    # to ensure the correlation function have the same definition of correlation
+    # lenght as is used in RMS and gaussianfft algorithm.
+    print(f"Variogram name:  {variogram_name}")
+    if variogram_name.upper() == "GAUSSIAN":
+        model = gs.Gaussian(
+            dim=3,
+            var=1.0,
+            len_scale=[xrange, yrange, zrange],
+            angles=np.pi * (0.5 - azimuth / 180.0),
+            rescale=math.sqrt(3),
+        )
+    elif variogram_name.upper() == "EXPONENTIAL":
+        model = gs.Exponential(
+            dim=3,
+            var=1.0,
+            len_scale=[xrange, yrange, zrange],
+            angles=np.pi * (0.5 - azimuth / 180.0),
+            rescale=3,
+        )
+    elif variogram_name.upper() == "GENERAL_EXPONENTIAL":
+        model = gs.Stable(
+            dim=3,
+            var=1.0,
+            len_scale=[xrange, yrange, zrange],
+            angles=np.pi * (0.5 - azimuth / 180.0),
+            rescale=math.pow(3.0, 1 / alpha),
+        )
+    else:
+        raise ValueError(f"Unknown variogram type: {variogram_name} ")
+    print(f"Start seed: {start_seed}")
+
+    print(f"Simulate field with size: nx={nx},ny={ny} nz={nz} ")
+    srf = gs.SRF(model, seed=start_seed)
+    field_srf = srf.structured([x, y, z], store="Field")
+    # print(f"Field: {srf.field_names} ")
+    # print(f"Field shape: {srf.field_shape} ")
+    # 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")
+    if settings["grid_size"]["use_eclipse_grid_index_origo"]:
+        field_result = np.zeros((nx, ny, nz), dtype=np.float32)
+        j_indices = -np.arange(ny) + ny - 1
+        field_result[:, j_indices, :] = field[:, :, :]
+        return field_result
+
+    return field
+
+
 def export_field(settings, field3D):
     # Export initial ensemble field
     nx, ny, nz = settings["field"]["grid_dimension"]

tests/jobs/localisation/example_case/sim_field_local.ert

@@ -2,6 +2,7 @@ DEFINE <USER>              $USER
 DEFINE <SCRATCH>           /scratch/fmu
 DEFINE <CASE_DIR>          sim_field_local
 DEFINE <ENSEMBLE_SEED_FILE> randomseeds.txt
+DEFINE <TEST_SETTINGS>  <CONFIG_PATH>/example_case.yml
 INSTALL_JOB  SIM_FIELD  scripts/FM_SIM_FIELD
 -----------------------------------------------------
 -- Observations