Casegen_UpCaRs: Add option to shift origin position to anywhere in the model

src/subscript/casegen_upcars/casegen_upcars.py

@@ -161,7 +161,10 @@ def main():
         centroid_x = get_value(geometry["centroid_x"], parser.centroid_x)
         centroid_y = get_value(geometry["centroid_y"], parser.centroid_y)
         origin_x = get_value(geometry.get("origin_x", 0.0), parser.origin_x)
+        origin_x_pos = get_value(geometry.get("origin_x_pos", 0.0), parser.origin_x_pos)
         origin_y = get_value(geometry.get("origin_y", 0.0), parser.origin_y)
+        origin_y_pos = get_value(geometry.get("origin_y_pos", 0.0), parser.origin_y_pos)
+        origin_top = get_value(geometry.get("origin_top", 0.0), parser.origin_top)
         rotation = get_value(geometry.get("rotation", 0.0), parser.rotation)
 
         # Merge streak intro background matrix
@@ -419,6 +422,9 @@ def main():
         origin_x,
         origin_y,
         rotation,
+        origin_x_pos,
+        origin_y_pos,
+        origin_top,
         fracture_length_x,
         fracture_offset_x,
         fracture_height_x,

src/subscript/casegen_upcars/model.py

@@ -62,6 +62,9 @@ def __init__(
         origin_x=0.0,
         origin_y=0.0,
         rotation=0.0,
+        origin_x_pos=0.0,
+        origin_y_pos=0.0,
+        origin_top=0.0,
         fracture_length_x=1.0,
         fracture_offset_x=0.0,
         fracture_height_x=1.0,
@@ -156,6 +159,9 @@ def __init__(
         self._origin_x = origin_x
         self._origin_y = origin_y
         self._rotation = rotation
+        self._origin_x_pos = origin_x_pos
+        self._origin_y_pos = origin_y_pos
+        self._origin_top = origin_top
 
         self._n_faults_x = self._matrix_x_count + (1 if fracture_at_boundary else -1)
         self._n_faults_y = self._matrix_y_count + (1 if fracture_at_boundary else -1)
@@ -337,34 +343,37 @@ def __init__(
 
         self._throws = []
 
-        print("Initializing model")
-        print(f"  LX x LY x LZ: {self._lx} x {self._ly} x {self._lz}")
-        print(f"  NX x NY x NZ: {self._total_nx} x {self._total_ny} x {self._total_nz}")
-        print(f"  dx x dy x dz: {dx} x {dy} x {dz}")
-
-        print(f"  Matrix Element in X-direction: {nMatrixY}")
-        print(f"  Matrix Element in Y-direction: {nMatrixX}")
-        print(f"  Layers: {nz}")
-        print(f"  Shape Factor: {radius_x} (a), {radius_y} (b), {radius_z} (c)")
-        print(f"  Top geometry: {top}")
-        print(f"  Model origin: {origin_x}, {origin_y}")
-        print(f"  Model coordinate rotation: {rotation}°")
-        print(f"  Tilting: {tilt}")
-        print("  Fracture:")
-        print(f"    Thickness: {fractureThickness}")
-        print(f"    Cell Count: {fracture_cell_count}")
-        print(f"    At Boundary: {fracture_at_boundary}")
-        print("    X-dir fractures:")
-        print(f"      Length: {fracture_length_x}")
-        print(f"      Offset: {fracture_offset_x}")
-        print(f"      Height: {fracture_height_x}")
-        print(f"      Z-offset: {fracture_zoffset_x}")
-        print("    Y-dir fractures:")
-        print(f"      Length: {fracture_length_y}")
-        print(f"      Offset: {fracture_offset_y}")
-        print(f"      Height: {fracture_height_y}")
-        print("      Z-offset: {fracture_zoffset_y}")
-
+        print(
+            f"""
+Initializing model
+  LX x LY x LZ: {self._lx} x {self._ly} x {self._lz}
+  NX x NY x NZ: {self._total_nx} x {self._total_ny} x {self._total_nz}
+  dx x dy x dz: {dx} x {dy} x {dz}
+  Matrix Element in X-direction: {nMatrixY}
+  Matrix Element in Y-direction: {nMatrixX}
+  Layers: {nz}
+  Shape Factor: {radius_x} (a), {radius_y} (b), {radius_z} (c)
+  Top geometry: {top}
+  Model origin top: {origin_top}
+  Model origin: {origin_x}, {origin_y} at {origin_x_pos}, {origin_y_pos}
+  Model coordinate rotation: {rotation}°
+  Tilting: {tilt}
+  Fracture:
+    Thickness: {fractureThickness}
+    Cell Count: {fracture_cell_count}
+    At Boundary: {fracture_at_boundary}
+    X-dir fractures:
+      Length: {fracture_length_x}
+      Offset: {fracture_offset_x}
+      Height: {fracture_height_x}
+      Z-offset: {fracture_zoffset_x}
+    Y-dir fractures:
+      Length: {fracture_length_y}
+      Offset: {fracture_offset_y}
+      Height: {fracture_height_y}
+      Z-offset: {fracture_zoffset_y}
+"""
+        )
         self.dict_info["nx"] = self._total_nx
         self.dict_info["ny"] = self._total_ny
         self.dict_info["nz"] = self._total_nz
@@ -468,6 +477,18 @@ def _build_grid(self):
         x_mid = self._centroid_x * self._lx
         y_mid = self._centroid_y * self._ly
 
+        rotation = np.radians(self._rotation)
+
+        origin_x_base = self._origin_x_pos * self._lx
+        origin_y_base = self._origin_y_pos * self._ly
+
+        origin_x_turn = (
+            math.cos(rotation) * origin_x_base + math.sin(rotation) * origin_y_base
+        )
+        origin_y_turn = (
+            -math.sin(rotation) * origin_x_base + math.cos(rotation) * origin_y_base
+        )
+
         cell_dx = np.full((self._total_nx, self._total_ny, self._total_nz), self._dx)
         cell_dy = np.full((self._total_nx, self._total_ny, self._total_nz), self._dy)
         cell_dz = np.empty((self._total_nx, self._total_ny, self._total_nz))
@@ -490,6 +511,15 @@ def _build_grid(self):
 
         self._xv, self._yv = np.meshgrid(self._x, self._y)
         if self._a * self._b * self._c != 0.0:
+            origin_z = -self._c * np.sqrt(
+                np.clip(
+                    1.0
+                    - (origin_x_base - x_mid) ** 2 / self._a**2
+                    - (origin_y_base - y_mid) ** 2 / self._b**2,
+                    0,
+                    None,
+                )
+            ) + (origin_x_base - x_mid) * math.tan(math.radians(self._tilt))
             self._zv = -self._c * np.sqrt(
                 np.clip(
                     1.0
@@ -500,21 +530,26 @@ def _build_grid(self):
                 )
             ) + (self._xv - x_mid) * math.tan(math.radians(self._tilt))
         else:
+            origin_z = (origin_x_base - x_mid) * math.tan(math.radians(self._tilt))
             self._zv = (self._xv - x_mid) * math.tan(math.radians(self._tilt))
 
-        rotation = np.radians(self._rotation)
-        rotation_matrix = np.array(
+        rotationMatrix = np.array(
             [
                 [np.cos(rotation), np.sin(rotation)],
                 [-np.sin(rotation), np.cos(rotation)],
             ]
         )
         self._xv, self._yv = np.einsum(
-            "ji, mni -> jmn", rotation_matrix, np.dstack([self._xv, self._yv])
+            "ji, mni -> jmn", rotationMatrix, np.dstack([self._xv, self._yv])
         )
-        self._xv += self._origin_x
-        self._yv += self._origin_y
-        self._zv += self._top - self._zv.min()
+
+        self._xv += self._origin_x - origin_x_turn
+        self._yv += self._origin_y - origin_y_turn
+
+        if self._origin_top > 0:
+            self._zv += self._origin_top - origin_z
+        else:
+            self._zv += self._top - self._zv.min()
 
     def distribute_property(self):
         """

src/subscript/casegen_upcars/udf_arg_parser.py

@@ -569,6 +569,15 @@ def fill_parser(parser):
         help="Origin coordinate of model in X-direction",
     )
 
+    parser.add_argument(
+        "--originXPos",
+        "--origin_x_pos",
+        type=float,
+        dest="origin_x_pos",
+        required=False,
+        help="Origin position as fraction of model size in X-direction",
+    )
+
     parser.add_argument(
         "--originY",
         "--origin_y",
@@ -578,6 +587,24 @@ def fill_parser(parser):
         help="Origin coordinate of model in Y-direction",
     )
 
+    parser.add_argument(
+        "--originYPos",
+        "--origin_y_pos",
+        type=float,
+        dest="origin_y_pos",
+        required=False,
+        help="Origin position as fraction of model size in Y-direction",
+    )
+
+    parser.add_argument(
+        "--originTop",
+        "--origin_top",
+        type=float,
+        dest="origin_top",
+        required=False,
+        help="Origin top depth",
+    )
+
     parser.add_argument(
         "--rotation",
         type=float,

tests/test_casegen_upcars.py

@@ -235,6 +235,56 @@ def test_demo_large_scale_with_coordinate_transformation(tmp_path, mocker):
     assert data_frame.Values["rotation"] == 15.0
 
 
+def test_demo_large_scale_with_origin_shifting(tmp_path, mocker):
+    """Test casegen_upcars on demo_large_scale.yaml with coordinate transformation"""
+    shutil.copytree(DATADIR, tmp_path / TESTDATA)
+    os.chdir(tmp_path / TESTDATA)
+
+    base_name = "TEST_LARGE_WITH_ORIGIN_SHIFTING"
+    mocker.patch(
+        "sys.argv",
+        [
+            "casegen_upcars",
+            "demo_large_scale.yaml",
+            "--et",
+            "dump_value.tmpl",
+            "--origin_x_pos",
+            "0.1",
+            "--origin_y_pos",
+            "0.8",
+            "--origin_top",
+            "1000.0",
+            "--base",
+            base_name,
+        ],
+    )
+    casegen_upcars.main()
+
+    # check that all output files are generated
+    for pre, suf in zip(
+        ["", "fipnum_", "gridinc_", "satnum_", "swat_"],
+        [".DATA", ".INC", ".GRDECL", ".INC", ".INC"],
+    ):
+        assert Path(pre + base_name + suf).exists()
+        if suf != ".DATA":
+            assert opm.io.Parser().parse(str(pre + base_name + suf))
+
+    # check some key parameters in output file
+    data_frame = pd.read_csv(base_name + ".DATA", index_col=0)
+    assert data_frame.Values["nx"] == 77
+    assert data_frame.Values["ny"] == 72
+    assert data_frame.Values["nz"] == 27
+    assert data_frame.Values["lx"] == 7700.0
+    assert data_frame.Values["ly"] == 7200.0
+    assert data_frame.Values["lz"] == 355.0
+    assert data_frame.Values["poro"] == 0.1711
+    assert data_frame.Values["originX"] == 0.0
+    assert data_frame.Values["originY"] == 0.0
+    assert data_frame.Values["rotation"] == 0.0
+    assert data_frame.Values["top"] == 1000.0
+    assert data_frame.Values["bottom"] == 1355.0
+
+
 def test_demo_large_scale_with_cmdline_streaks(tmp_path, mocker):
     """Test casegen_upcars on demo_large_scale.yaml with some streaks"""
     shutil.copytree(DATADIR, tmp_path / TESTDATA)

tests/testdata_casegen_upcars/demo_large_scale.yaml

@@ -48,7 +48,12 @@ Geometry:
   origin_y: 0
   # Coordinate rotation angle in degree ('--rotation')
   rotation: 0
-
+  # The location of origin as fraction of the model lx (0 = West most - Default, 1 = East most) ('--origin_x_pos', '--originXPos')
+  origin_x_pos: 1.0
+  # The location of origin as fraction of the model ly (0 = South most - Default, 1 = North most) ('--origin_y_pos', '--originYPos')
+  origin_y_pos: 1.0
+  # The top of origin point, default = 0 in which case it's ignored and model will use top keyword, if > 0, CaseGen will ensure that the origin point will start at this depth ('--origin_top', '--originTop')
+  origin_top: 0.0
 Layers:
   Background Matrix:
     NZ : 27                   # Total number of cells in z-direction, including heterogeneous layers ('--matrix_nz)

tests/testdata_casegen_upcars/demo_small_scale.yaml

@@ -48,7 +48,12 @@ Geometry:
   origin_y: 0
   # Coordinate rotation angle in degree ('--rotation')
   rotation: 0
-
+  # The location of origin as fraction of the model lx (0 = West most - Default, 1 = East most) ('--origin_x_pos', '--originXPos')
+  origin_x_pos: 1.0
+  # The location of origin as fraction of the model ly (0 = South most - Default, 1 = North most) ('--origin_y_pos', '--originYPos')
+  origin_y_pos: 1.0
+  # The top of origin point, default = 0 in which case it's ignored and model will use top keyword, if > 0, CaseGen will ensure that the origin point will start at this depth ('--origin_top', '--originTop')
+  origin_top: 0.0
 
 Layers:
   Background Matrix: