From b5dccd66d510d18ea37883a9976d86e51f5a29c1 Mon Sep 17 00:00:00 2001 From: David Landa Marban Date: Wed, 17 Jun 2026 14:26:27 +0200 Subject: [PATCH] Bug fixes for dual porosity --- .github/workflows/ci_pycopm_ubuntu.yml | 4 +- examples/decks/MODEL7.DATA | 111 ++++++++ examples/decks/MODEL8.DATA | 120 +++++++++ src/pycopm/core/pycopm.py | 4 +- src/pycopm/utils/generate_files.py | 17 +- src/pycopm/utils/mapping_methods.py | 127 ++++++--- tests/test_7_transmissibilities.py | 6 +- tests/test_8_dual.py | 339 +++++++++++++++++++++---- 8 files changed, 632 insertions(+), 96 deletions(-) create mode 100644 examples/decks/MODEL7.DATA create mode 100644 examples/decks/MODEL8.DATA diff --git a/.github/workflows/ci_pycopm_ubuntu.yml b/.github/workflows/ci_pycopm_ubuntu.yml index d6e5e17..1e8bece 100644 --- a/.github/workflows/ci_pycopm_ubuntu.yml +++ b/.github/workflows/ci_pycopm_ubuntu.yml @@ -49,13 +49,13 @@ jobs: - name: Check code style and linting run: | - black --target-version py312 --check src/ tests/ + black --target-version py312 src/ tests/ --check pylint src/ tests/ mypy --ignore-missing-imports src/ tests/ - name: Run tests run: | - pytest --cov=pycopm --exitfirst --cov-report term-missing --basetemp=test_outputs tests/ + pytest --cov=pycopm --cov-report term-missing --basetemp=test_outputs tests/ --exitfirst - name: Run hello world example run: | diff --git a/examples/decks/MODEL7.DATA b/examples/decks/MODEL7.DATA new file mode 100644 index 0000000..a70ce30 --- /dev/null +++ b/examples/decks/MODEL7.DATA @@ -0,0 +1,111 @@ +-- Copyright (C) 2026 NORCE Research AS +-- flow MODEL7.DATA +-- pycopm -i MODEL7.DATA -w STANDARD -t 2 -a max -z 1:2 -l S +-- pycopm -i MODEL7.DATA -w DUAL -t 2 -a max -z 1:2 -l D -dual 'poro == 0.1' +---------------------------------------------------------------------------- +RUNSPEC +---------------------------------------------------------------------------- +DIMENS +3 3 2 / + +EQLDIMS +/ + +TABDIMS +/ + +CO2STORE +GAS +WATER + +METRIC + +START +1 JAN 2025 / + +WELLDIMS +/ + +UNIFOUT +---------------------------------------------------------------------------- +GRID +---------------------------------------------------------------------------- +INIT + +DX +18*1E2 / + +DY +18*1E2 / + +DZ +18*1E2 / + +TOPS +9*0 / + +FLUXNUM +1 1 2 1 2 2 1 1 2 +2 1 1 2 2 1 2 1 1 +/ + +EQUALREG +PORO 0 0 F / +PORO 0.25 1 F / +PORO 1e-1 2 F / +PERMX 0 0 F / +PERMX 1.172745396974317 1 F / +PERMX 1.172745396974317e-2 2 F / +PERMY 0 0 F / +PERMY 1.172745396974317 1 F / +PERMY 1.172745396974317e-2 2 F / +PERMZ 0 0 F / +PERMZ 1.172745396974317 1 F / +PERMZ 1.172745396974317e-2 2 F / +/ + +NNC +1 1 1 2 1 2 3 / +1 1 1 1 2 2 7 / +2 2 2 3 2 1 11 / +/ +---------------------------------------------------------------------------- +PROPS +---------------------------------------------------------------------------- +SGWFN +0 0 1 0 +1 1 0 0 / + +ROCK +276 4.934e-05 / +---------------------------------------------------------------------------- +SOLUTION +---------------------------------------------------------------------------- +EQUIL +0 200 0 0 0 0 1 1 0 / + +RPTRST +BASIC=2 RPORV / +---------------------------------------------------------------------------- +SUMMARY +---------------------------------------------------------------------------- +BPR +3 3 2 / +/ +---------------------------------------------------------------------------- +SCHEDULE +---------------------------------------------------------------------------- +RPTRST +BASIC=2 RPORV / + +SOURCE +1 1 1 WATER 1E+04 / +/ +TSTEP +9131.25 / + +SOURCE +1 1 1 WATER 0 / +/ +TSTEP +3652.5 / diff --git a/examples/decks/MODEL8.DATA b/examples/decks/MODEL8.DATA new file mode 100644 index 0000000..55bbbd5 --- /dev/null +++ b/examples/decks/MODEL8.DATA @@ -0,0 +1,120 @@ +-- Copyright (C) 2026 NORCE Research AS +-- flow MODEL8.DATA +-- pycopm -i MODEL8.DATA -w STANDARD -t 2 -a max -z 1:2 -l S +-- pycopm -i MODEL8.DATA -w DUAL -t 2 -a max -z 1:2 -l D -dual 'poro == 0.1' +---------------------------------------------------------------------------- +RUNSPEC +---------------------------------------------------------------------------- +DIMENS +3 2 2 / + +EQLDIMS +/ + +TABDIMS +/ + +CO2STORE +GAS +WATER + +METRIC + +START +1 JAN 2025 / + +WELLDIMS +/ + +UNIFOUT +---------------------------------------------------------------------------- +GRID +---------------------------------------------------------------------------- +INIT + +COORD + 0 0 0 0 0 2E3 + 2.5E3 0 0 2.5E3 0 2E3 + 5E3 0 0 5E3 0 2E3 + 7.5E3 0 0 7.5E3 0 2E3 + 0 1E3 0 0 1E3 2E3 + 2.5E3 1E3 0 2.5E3 1E3 2E3 + 5E3 1E3 0 5E3 1E3 2E3 + 7.5E3 1E3 0 7.5E3 1E3 2E3 + 0 2E3 0 0 2E3 2E3 + 2.5E3 2E3 0 2.5E3 2E3 2E3 + 5E3 2E3 0 5E3 2E3 2E3 + 7.5E3 2E3 0 7.5E3 2E3 2E3 / + +ZCORN + 5E2 5E2 0 0 -5E2 -5E2 5E2 5E2 0 0 -5E2 -5E2 + 5E2 5E2 -2.5E2 -2.5E2 -5E2 -5E2 5E2 5E2 -2.5E2 -2.5E2 -5E2 -5E2 + +15E2 15E2 7.5E2 7.5E2 2.5E2 2.5E2 15E2 15E2 7.5E2 7.5E2 2.5E2 2.5E2 +15E2 15E2 5E2 5E2 2.5E2 2.5E2 15E2 15E2 5E2 5E2 2.5E2 2.5E2 + +15E2 15E2 7.5E2 7.5E2 2.5E2 2.5E2 15E2 15E2 7.5E2 7.5E2 2.5E2 2.5E2 +15E2 15E2 5E2 5E2 2.5E2 2.5E2 15E2 15E2 5E2 5E2 2.5E2 2.5E2 +15E2 15E2 10E2 10E2 10E2 10E2 15E2 15E2 10E2 10E2 10E2 10E2 +15E2 15E2 10E2 10E2 10E2 10E2 15E2 15E2 10E2 10E2 10E2 10E2 / + +FLUXNUM +1 2 2 0 2 0 +1 1 1 0 2 0 +/ + +EQUALREG +PORO 0 0 F / +PORO 0.25 1 F / +PORO 1e-1 2 F / +PERMX 0 0 F / +PERMX 1.172745396974317 1 F / +PERMX 1.172745396974317e-2 2 F / +PERMY 0 0 F / +PERMY 1.172745396974317 1 F / +PERMY 1.172745396974317e-2 2 F / +PERMZ 0 0 F / +PERMZ 1.172745396974317 1 F / +PERMZ 1.172745396974317e-2 2 F / +/ + +---------------------------------------------------------------------------- +PROPS +---------------------------------------------------------------------------- +SGWFN +0 0 1 0 +1 1 0 0 / + +ROCK +276 4.934e-05 / +---------------------------------------------------------------------------- +SOLUTION +---------------------------------------------------------------------------- +EQUIL +-5E2 200 0 0 0 0 1 1 0 / + +RPTRST +BASIC=2 RPORV / +---------------------------------------------------------------------------- +SUMMARY +---------------------------------------------------------------------------- +BPR +3 1 2 / +/ +---------------------------------------------------------------------------- +SCHEDULE +---------------------------------------------------------------------------- +RPTRST +BASIC=2 RPORV / + +SOURCE +1 1 1 WATER 1E+04 / +/ +TSTEP +9131.25 / + +SOURCE +1 1 1 WATER 0 / +/ +TSTEP +3652.5 / diff --git a/src/pycopm/core/pycopm.py b/src/pycopm/core/pycopm.py index 03a4cbb..eefe860 100644 --- a/src/pycopm/core/pycopm.py +++ b/src/pycopm/core/pycopm.py @@ -374,7 +374,9 @@ def load_parser(): "--dual", default="", help="Set the criterium to differentiate net and non-net in coarsening using a static " - "variable, e.g., 'poro <= 0.1' ('' by default).", + "variable. To remove the vertical transfer function (FT) between net and not-net cells, " + "add to the command ', vertical TF = 0', e.g., 'poro <= 0.1' (which includes vertical TF) " + "or 'poro <= 0.1, vertical TF = 0' ('' by default)", ) return vars(parser.parse_known_args()[0]) diff --git a/src/pycopm/utils/generate_files.py b/src/pycopm/utils/generate_files.py index 1385465..dc653bf 100644 --- a/src/pycopm/utils/generate_files.py +++ b/src/pycopm/utils/generate_files.py @@ -670,22 +670,25 @@ def initialize_variables(dic): dic["facpermz"] = dic["ini"]["PERMZ"][i] / float(val) break dic["msk"] = np.ones(len(dic["porv"])) + dic["dual_vertical_tf"] = True if dic["dual"]: dual = dic["dual"].split(" ") + val = dual[2].replace(",", "") + dic["dual_vertical_tf"] = "vertical TF = 0" not in dic["dual"] quan = np.array(dic["ini"][dual[0].upper()]) dic["dual"] = 1 if dual[1] == "==": - dic["msk"][dic["porv"] > 0] = quan != float(dual[2]) + dic["msk"][dic["porv"] > 0] = quan != float(val) elif dual[1] == ">=": - dic["msk"][dic["porv"] > 0] = quan < float(dual[2]) + dic["msk"][dic["porv"] > 0] = quan < float(val) elif dual[1] == "<=": - dic["msk"][dic["porv"] > 0] = quan > float(dual[2]) + dic["msk"][dic["porv"] > 0] = quan > float(val) elif dual[1] == "<": - dic["msk"][dic["porv"] > 0] = quan >= float(dual[2]) + dic["msk"][dic["porv"] > 0] = quan >= float(val) elif dual[1] == ">": - dic["msk"][dic["porv"] > 0] = quan <= float(dual[2]) + dic["msk"][dic["porv"] > 0] = quan <= float(val) elif dual[1] == "!=": - dic["msk"][dic["porv"] > 0] = quan == float(dual[2]) + dic["msk"][dic["porv"] > 0] = quan == float(val) else: print(f"Unknow criterium for non-net cells ({dic['dual']}).") sys.exit() @@ -915,7 +918,7 @@ def write_props(dic, n): tmp = compact_format(dic[f"{name}_c"]) if ( "*" in tmp[0] - and not (dic["trans"] > 0 and name in ["tranx", "trany"]) + and not (dic["trans"] > 0 and name in ["tranx", "trany", "tranz"]) and (not dic["dual"] or n > dic["ntot"]) ): if int(tmp[0].split("*")[0]) == n: diff --git a/src/pycopm/utils/mapping_methods.py b/src/pycopm/utils/mapping_methods.py index cf9cca4..d06dc99 100644 --- a/src/pycopm/utils/mapping_methods.py +++ b/src/pycopm/utils/mapping_methods.py @@ -146,12 +146,38 @@ def map_properties(dic, actnum, z_t, z_b, z_b_t, v_c): for name in dic["props"] + dic["rptrst"]: bar_animation() q = dic[name].copy() - q[dic["msk"] == 0] = 0.0 - q1 = dic[name].copy() - q1[dic["msk"] == 1] = 0.0 + if dic["dual"]: + q1 = dic[name].copy() + q[dic["msk"] == 0] = 0.0 + q1[dic["msk"] == 1] = 0.0 + if name == "tranx": + msk = np.roll( + dic["msk"].reshape((dic["zn"], dic["yn"], dic["xn"])), + -1, + axis=2, + ).ravel() + q[msk == 0] = 0.0 + q1[msk == 1] = 0.0 + elif name == "trany": + msk = np.roll( + dic["msk"].reshape((dic["zn"], dic["yn"], dic["xn"])), + -1, + axis=1, + ).ravel() + q[msk == 0] = 0.0 + q1[msk == 1] = 0.0 + elif name == "tranz": + msk = np.roll( + dic["msk"].reshape((dic["zn"], dic["yn"], dic["xn"])), + -1, + axis=0, + ).ravel() + q[msk == 0] = 0.0 + q1[msk == 1] = 0.0 if not dic["show"]: c_c = pd.Series(q).groupby(dic["con"]).sum() - c_c1 = pd.Series(q1).groupby(dic["con"]).sum() + if dic["dual"]: + c_c1 = pd.Series(q1).groupby(dic["con"]).sum() if name in ["permx", "permy"]: dic[f"{name}_c"] = [ f"{val/h_t:E}" if h_t * val > 0 else "0" @@ -193,17 +219,11 @@ def map_properties(dic, actnum, z_t, z_b, z_b_t, v_c): ] else: dic[f"{name}_c"] = [ - f"{val*l_a/l_t:E}" if val > 0 else "0" - for val, l_a, l_t in zip( - c_c, dic[f"{drc}_a"], dic[f"{drc}_tot"] - ) + f"{val:E}" if val > 0 else "0" for val in c_c ] if dic["dual"]: dic[f"{name}_dual_c"] = [ - f"{val*l_a/l_t:E}" if val > 0 else "0" - for val, l_a, l_t in zip( - c_c1, dic[f"{drc}_a"], dic[f"{drc}_tot"] - ) + f"{val:E}" if val > 0 else "0" for val in c_c1 ] elif name == "trany": if "y" in dic["coardir"]: @@ -232,17 +252,11 @@ def map_properties(dic, actnum, z_t, z_b, z_b_t, v_c): ] else: dic[f"{name}_c"] = [ - f"{val*l_a/l_t:E}" if val > 0 else "0" - for val, l_a, l_t in zip( - c_c, dic[f"{drc}_a"], dic[f"{drc}_tot"] - ) + f"{val:E}" if val > 0 else "0" for val in c_c ] if dic["dual"]: dic[f"{name}_dual_c"] = [ - f"{val*l_a/l_t:E}" if val > 0 else "0" - for val, l_a, l_t in zip( - c_c1, dic[f"{drc}_a"], dic[f"{drc}_tot"] - ) + f"{val:E}" if val > 0 else "0" for val in c_c1 ] elif name == "tranz": if "z" in dic["coardir"]: @@ -263,11 +277,28 @@ def map_properties(dic, actnum, z_t, z_b, z_b_t, v_c): ] else: c_ls = pd.Series(q).groupby(dic["con"]).last() - d_ls = pd.Series(dic["d_z"]).groupby(dic["con"]).last() + d_ls = pd.Series(v_c).groupby(dic["con"]).last() + d_ld = pd.Series(v_c).groupby(dic["con"]).first() + d_ld = list(d_ld[dic["xn"] * dic["yn"] :]) + list( + d_ld[: dic["xn"] * dic["yn"]] + ) + l_dd = list(v_tot[dic["xn"] * dic["yn"] :]) + list( + v_tot[: dic["xn"] * dic["yn"]] + ) dic[f"{name}_c"] = [ - f"{val*d_l/l_t:E}" if val > 0 else "0" - for val, d_l, l_t in zip(c_ls, d_ls, dic[f"{drc}_tot"]) + f"{val*(d_l+d_d)/(l_t + l_d):E}" if val > 0 else "0" + for val, d_l, d_d, l_t, l_d in zip( + c_ls, d_ls, d_ld, v_tot, l_dd + ) ] + if dic["dual"]: + c_ls = pd.Series(q1).groupby(dic["con"]).last() + dic[f"{name}_dual_c"] = [ + f"{val*(d_l+d_d)/(l_t + l_d):E}" if val > 0 else "0" + for val, d_l, d_d, l_t, l_d in zip( + c_ls, d_ls, d_ld, v_tot, l_dd + ) + ] else: dic[f"{name}_c"] = [ f"{val*l_a/l_t:E}" if val > 0 else "0" @@ -1384,16 +1415,7 @@ def handle_cp_grid(dic): def handle_dual(dic): - """ - Create the dual grid by extending it in the j-direction - - Args: - dic (dict): Global dictionary - - Returns: - dic (dict): Modified global dictionary - - """ + """Create the dual grid by extending it in the j-direction""" nxy = dic["nx"] * dic["ny"] dic["cr"] = list(dic["cr"]) dy = 1.075 * max( @@ -1426,6 +1448,7 @@ def handle_dual(dic): and float(dic["porv_dual_c"][idd]) > 0 and float(dic["porv_c"][idd]) > 0 and not np.isnan(dic["tranz"][ind]) + and dic["dual_vertical_tf"] ): trnz += dic["tranz"][ind] if i < dic["xn"] - 1: @@ -1461,6 +1484,40 @@ def handle_dual(dic): ): trnmy += dic["trany"][ind - dic["xn"]] k += 1 + ind = i + j * dic["xn"] + k * nxy + idd = i + j * dic["xn"] + n * nxy + if i < dic["xn"] - 1: + if ( + dic["msk"][ind] == 1 + and dic["msk"][ind + 1] == 0 + and float(dic["porv_dual_c"][idd + 1]) > 0 + and not np.isnan(dic["tranx"][ind]) + ): + trnx += dic["tranx"][ind] + if 0 < i: + if ( + dic["msk"][ind - 1] == 0 + and dic["msk"][ind] == 1 + and float(dic["porv_dual_c"][idd - 1]) > 0 + and not np.isnan(dic["tranx"][ind - 1]) + ): + trnmx += dic["tranx"][ind - 1] + if j < dic["yn"] - 1: + if ( + dic["msk"][ind] == 1 + and dic["msk"][ind + dic["xn"]] == 0 + and float(dic["porv_dual_c"][idd + dic["xn"]]) > 0 + and not np.isnan(dic["trany"][ind]) + ): + trny += dic["trany"][ind] + if 0 < j: + if ( + dic["msk"][ind - dic["xn"]] == 0 + and dic["msk"][ind] == 1 + and float(dic["porv_dual_c"][idd - dic["xn"]]) > 0 + and not np.isnan(dic["trany"][ind - dic["xn"]]) + ): + trnmy += dic["trany"][ind - dic["xn"]] k += 1 if trnz > 0: dic[ @@ -1509,12 +1566,10 @@ def handle_dual(dic): tmp += [f"{default}" for _ in range(dic["nx"])] if prop == "fluxnum": tmp += tmp2 - elif prop in ["porv", "poro", "tranx", "trany"]: + elif prop in ["porv", "poro", "tranx", "trany", "tranz"]: tmp += dic[f"{prop}_dual_c"][i * nxy : (i + 1) * nxy] elif prop in ["permx", "permy", "permz"]: tmp += tmp0 - elif prop in ["tranz"]: - tmp += [f"{default}" for _ in range(nxy)] else: tmp += dic[f"{prop}_c"][i * nxy : (i + 1) * nxy] dic[f"{prop}_c"] = tmp diff --git a/tests/test_7_transmissibilities.py b/tests/test_7_transmissibilities.py index b7eee15..29d7aa6 100644 --- a/tests/test_7_transmissibilities.py +++ b/tests/test_7_transmissibilities.py @@ -15,8 +15,8 @@ def test_7_transmissibilities(flow, tmp_path, monkeypatch): repo_root = Path(__file__).parents[1] monkeypatch.chdir(tmp_path) values = [ - [696.14886, 774.826, 5008.1323, 13739.814453], - [696.14886, 774.826, 0, 13739.814453], + [827.3911743164062, 968.17822265625, 5008.1323, 13739.814453], + [827.3911743164062, 968.17822265625, 8564.0341796875, 13739.814453], ] for i, sub in enumerate(["1", "2"]): subprocess.run( @@ -45,7 +45,7 @@ def test_7_transmissibilities(flow, tmp_path, monkeypatch): for j, n in enumerate(["X", "Y", "Z", "NNC"]): assert ( abs(np.sum(np.array(init[f"TRAN{n}"])) - values[i][j]) < 1e-5 - ), f"Issue in TRAN{n} with -t {i}" + ), f"Issue in TRAN{n} with -t {sub}" subprocess.run( [ "pycopm", diff --git a/tests/test_8_dual.py b/tests/test_8_dual.py index 2ef16bc..5cddd1b 100644 --- a/tests/test_8_dual.py +++ b/tests/test_8_dual.py @@ -1,63 +1,308 @@ # SPDX-FileCopyrightText: 2026 NORCE Research AS # SPDX-License-Identifier: GPL-3.0 +# pylint: disable=R0913,R0917 """Test the dual coarsening, i.e., splitting net and non-net cells""" from pathlib import Path import subprocess +import pytest import numpy as np from opm.io.ecl import EclFile as OpmFile -from opm.io.ecl import ERst as OpmRestart +TOL = {"rtol": 1e-5, "atol": 1e-8} -def test_8_dual(flow, tmp_path, monkeypatch): - """See examples/decks/MODEL6.DATA""" - repo_root = Path(__file__).parents[1] - monkeypatch.chdir(tmp_path) - subprocess.run( - [ - "pycopm", - "-i", - f"{repo_root}/examples/decks/MODEL6.DATA", - "-f", - flow, - "-z", + +@pytest.mark.parametrize( + "deck,z_range,dual_expr,expected", + [ + ( + "MODEL6.DATA", + "1:4", + "poro <= 0.1, vertical TF = 0", + { + "PORV": np.array( + [ + 6.25e8, + 1.25e9, + 6.25e8, + 6.25e8, + 6.25e8, + 6.25e8, + 1.25e9, + 6.25e8, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 2.50e8, + 2.50e8, + 2.50e8, + 2.50e8, + 0, + 0, + ] + ), + "TRANX": np.array( + [ + 3.456000e2, + 3.456000e2, + 3.456000e2, + 3.456000e2, + 3.456000e2, + 3.456000e2, + 3.456000e2, + 0, + 3.410807e-2, + 3.410807e-2, + 3.410807e-2, + 0, + ] + ), + "TRANY": np.zeros(12), + "TRANNNC": np.array([0.0682094, 0.0682094]), + "NNC1": np.array([2, 7]), + "NNC2": np.array([19, 22]), + }, + ), + ( + "MODEL6.DATA", "1:4", - "-w", - "DUAL", - "-dual", "poro <= 0.1", - "D", - "-t", - "2", - "-a", - "max", - ], - check=True, + { + "PORV": np.array( + [ + 6.25e8, + 1.25e9, + 6.25e8, + 6.25e8, + 6.25e8, + 6.25e8, + 1.25e9, + 6.25e8, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 2.50e8, + 2.50e8, + 2.50e8, + 2.50e8, + 0, + 0, + ] + ), + "TRANX": np.array( + [ + 3.456000e2, + 3.456000e2, + 3.456000e2, + 3.456000e2, + 3.456000e2, + 3.456000e2, + 3.456000e2, + 0, + 3.410807e-2, + 3.410807e-2, + 3.410807e-2, + 0, + ] + ), + "TRANY": np.zeros(12), + "TRANNNC": np.array( + [0.0682094, 2.5578527, 2.5578527, 2.5578527, 2.5578527, 0.0682094] + ), + "NNC1": np.array([2, 3, 4, 5, 6, 7]), + "NNC2": np.array([19, 19, 20, 21, 22, 22]), + }, + ), + ( + "MODEL7.DATA", + "1:2", + "poro == 0.1", + { + "PORV": np.array( + [ + 250000, + 500000, + 250000, + 250000, + 0, + 250000, + 250000, + 500000, + 250000, + 0, + 0, + 0, + 100000, + 0, + 100000, + 100000, + 200000, + 100000, + 100000, + 0, + 100000, + ] + ), + "TRANX": np.array( + [ + 4.000002e0, + 1.000002e0, + 0, + 0, + 0, + 1.000002e0, + 1.000002e0, + 0, + 0, + 0, + 1.000002e-2, + 1.101000e1, + 0, + 0, + 0, + ] + ), + "TRANY": np.array( + [ + 1.000002, + 0, + 1.000002, + 1.000002, + 1.000002, + 0, + 0, + 0, + 0.01000002, + 0.01000002, + 0.01000002, + 0, + 0.01000002, + 0, + 0, + ] + ), + "TRANNNC": np.array( + [ + 0.01980202, + 7, + 0.01980202, + 0.01980202, + 0.03960404, + 0.01980202, + 0.01980202, + 0.01980202, + 0.01980202, + 0.01980202, + 0.01980202, + 0.03960404, + 0.01980202, + 0.01980202, + 0.01980202, + ] + ), + "NNC1": np.array([1, 1, 2, 2, 2, 3, 4, 4, 6, 6, 7, 8, 8, 8, 9]), + "NNC2": np.array( + [13, 16, 13, 15, 17, 15, 16, 17, 17, 18, 19, 17, 19, 21, 21] + ), + }, + ), + ( + "MODEL8.DATA", + "1:2", + "poro == 0.1", + { + "PORV": np.array( + [ + 6.2500e8, + 1.5625e8, + 4.6875e8, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 1.8750e8, + 1.8750e8, + 0, + 3.1250e8, + 0, + ] + ), + "TRANX": np.array([1.000002, 1.000002, 0, 0.01000002, 0, 0.0]), + "TRANY": np.array([0, 0, 0, 0.18750042, 0, 0.0]), + "TRANNNC": np.array( + [0.01980202, 0.66445315, 0.12376262, 0.03960404, 0.6600673] + ), + "NNC1": np.array([1, 2, 2, 3, 3]), + "NNC2": np.array([11, 11, 14, 11, 12]), + }, + ), + ], +) +def test_8_dual(flow, tmp_path, monkeypatch, deck, z_range, dual_expr, expected): + """See examples/decks/MODEL7.DATA and MODEL.DATA""" + repo_root = Path(__file__).resolve().parents[1] + monkeypatch.chdir(tmp_path) + + cmd = [ + "pycopm", + "-i", + str(repo_root / f"examples/decks/{deck}"), + "-f", + flow, + "-z", + z_range, + "-w", + "DUAL", + "-dual", + dual_expr, + "D", + "-t", + "2", + "-m", + "all", + "-a", + "max", + ] + + subprocess.run(cmd, check=True) + + init = OpmFile("DUAL.INIT") + egrid = OpmFile("DUAL.EGRID") + + np.testing.assert_allclose( + init["PORV"], expected["PORV"], err_msg="PORV mismatch", **TOL + ) + np.testing.assert_allclose( + init["TRANX"], expected["TRANX"], err_msg="TRANX mismatch", **TOL ) - subprocess.run( - [ - flow, - "DUAL.DATA", - ], - check=True, + np.testing.assert_allclose( + init["TRANY"], expected["TRANY"], err_msg="TRANY mismatch", **TOL + ) + np.testing.assert_allclose( + init["TRANNNC"], expected["TRANNNC"], err_msg="TRANNNC mismatch", **TOL + ) + + np.testing.assert_array_equal( + egrid["NNC1"], expected["NNC1"], err_msg="NNC1 mismatch" ) - subprocess.run( - [ - flow, - f"{repo_root}/examples/decks/MODEL6.DATA", - f"--output-dir={tmp_path}", - ], - check=True, + np.testing.assert_array_equal( + egrid["NNC2"], expected["NNC2"], err_msg="NNC2 mismatch" ) - initr = OpmFile("MODEL6.INIT") - initd = OpmFile("DUAL.INIT") - assert ( - abs(np.sum(np.array(initr["PORV"])) - np.sum(np.array(initd["PORV"]))) < 1e-12 - ), "Issue with PORV and dual" - assert (len(initd["PORV"])) == 24, "Issue with dimensions of the dual model" - unrstr = OpmRestart("MODEL6.UNRST") - unrstd = OpmRestart("DUAL.UNRST") - assert ( - abs(unrstr["PRESSURE", 1][7,] - unrstd["PRESSURE", 1][7]) < 1e-1 - ), "Issue with the dual transsmissibility"