Merge pull request #78 from daavid00/dev

Bug fix for coarsening and -t 1 or -t 2 for models with TRANNNC

Author: daavid00

examples/decks/MODEL5.DATA

@@ -0,0 +1,95 @@
+-- This reservoir simulation deck is made available under the Open Database
+-- License: http://opendatacommons.org/licenses/odbl/1.0/. Any rights in
+-- individual contents of the database are licensed under the Database Contents
+-- License: http://opendatacommons.org/licenses/dbcl/1.0/
+-- Copyright (C) 2026 NORCE Research AS
+----------------------------------------------------------------------------
+RUNSPEC
+----------------------------------------------------------------------------
+DIMENS 
+4 1 2 /
+
+EQLDIMS
+/
+
+TABDIMS
+/
+
+CO2STORE
+GAS
+WATER
+
+METRIC
+
+START
+1 JAN 2025 /
+
+WELLDIMS
+/
+
+UNIFIN
+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
+  10E3   0 0   10E3   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
+  10E3 1E3 0   10E3 1E3 2E3 /
+
+ZCORN
+16*0
+10E2 10E2 4*5E2 10E2 10E2 10E2 10E2 4*5E2 10E2 10E2
+10E2 10E2 4*5E2 10E2 10E2 10E2 10E2 4*5E2 10E2 10E2
+16*10E2 /
+
+PORO
+8*1 /
+
+PERMX
+8*10.13 /
+
+COPY
+PERMX PERMY /
+PERMX PERMZ /
+/
+----------------------------------------------------------------------------
+PROPS
+----------------------------------------------------------------------------
+SGWFN
+0 0 1 0
+1 1 0 0 /
+
+ROCK
+276 4.934e-05 /
+----------------------------------------------------------------------------
+SOLUTION
+----------------------------------------------------------------------------
+EQUIL
+0 200.0 0 0 0 0 1 1 0 /
+
+RPTRST 
+BASIC=2 /
+----------------------------------------------------------------------------
+SUMMARY
+----------------------------------------------------------------------------
+----------------------------------------------------------------------------
+SCHEDULE
+----------------------------------------------------------------------------
+RPTRST
+BASIC=2 /
+
+SOURCE
+1 1 1 WATER 8.015165E+06 /
+/
+TSTEP
+1 /

src/pycopm/utils/generate_files.py

@@ -417,8 +417,10 @@ def create_deck(dic):
             print("\nThis is needed for the nnctrans, please wait.\n")
             os.chdir(dic["fol"])
             os.system(f"{dic['flow']} {dic['fol']}/{dic['write']}.DATA {dic['flags']}")
-            dic["ogrid"] = OpmFile(dic["write"] + ".EGRID")
-            if dic["ogrid"].count("NNC1") and dic["trans"] > 0:
+            if (
+                OpmFile(dic["write"] + ".EGRID").count("NNC1")
+                or OpmFile(dic["deck"] + ".EGRID").count("NNC1")
+            ) and dic["trans"] > 0:
                 handle_nnc_trans(dic)
             else:
                 print("\nNo nnctrans found.")
@@ -562,7 +564,10 @@ def initialize_variables(dic):
                     dic["mults"] += [name]
     for name in ["multnum", "fluxnum"]:
         if dic["ini"].count(name.upper()):
-            if np.max(dic["ini"][name.upper()]) > 1:
+            if (
+                np.max(dic["ini"][name.upper()]) > 1
+                or np.min(dic["ini"][name.upper()]) < 1
+            ):
                 dic["grids"] += [name]
     for name in ["thconr", "disperc"]:
         if dic["ini"].count(name.upper()):
@@ -633,51 +638,47 @@ def handle_nnc_trans(dic):
         dic (dict): Modified global dictionary
 
     """
-    nncdeck, withlines = [], False
-    with open(dic["write"] + ".DATA", "r", encoding=dic["encoding"]) as file:
-        for row in csv.reader(file):
-            nrwo = str(row)[2:-2].strip()
-            nncdeck.append(nrwo)
-            if withlines:
-                nncdeck.append("INCLUDE")
-                nncdeck.append(f"'{dic['label']}EDITNNC.INC' /\n")
-                withlines = False
-            elif nrwo == "EDIT":
-                if not dic["lines"]:
-                    nncdeck.append("INCLUDE")
-                    nncdeck.append(f"'{dic['label']}EDITNNC.INC' /\n")
-                else:
-                    withlines = True
-    with open(
-        f"{dic['fol']}/{dic['write']}.DATA",
-        "w",
-        encoding="utf8",
-    ) as file:
-        for row in nncdeck:
-            file.write(row + "\n")
     temp = OpmFile(dic["deck"] + ".EGRID")
     coa = OpmFile(dic["write"] + ".INIT")
-    coag = OpmFile(dic["write"] + ".EGRID")
+    editnnc = OpmFile(dic["write"] + ".EGRID").count("NNC1")
+    if editnnc:
+        nncdeck, indel, withlines = [], [], False
+        dic["coa_editnnc"] = []
+        dic["nnct"] = [[[] for _ in range(dic["yn"])] for _ in range(dic["xn"])]
+        with open(dic["write"] + ".DATA", "r", encoding=dic["encoding"]) as file:
+            for row in csv.reader(file):
+                nrwo = str(row)[2:-2].strip()
+                nncdeck.append(nrwo)
+                if withlines:
+                    nncdeck.append("INCLUDE")
+                    nncdeck.append(f"'{dic['label']}EDITNNC.INC' /\n")
+                    withlines = False
+                elif nrwo == "EDIT":
+                    if not dic["lines"]:
+                        nncdeck.append("INCLUDE")
+                        nncdeck.append(f"'{dic['label']}EDITNNC.INC' /\n")
+                    else:
+                        withlines = True
+        with open(
+            f"{dic['fol']}/{dic['write']}.DATA",
+            "w",
+            encoding="utf8",
+        ) as file:
+            for row in nncdeck:
+                file.write(row + "\n")
+        coag = OpmFile(dic["write"] + ".EGRID")
+        cnnc1 = np.array(coag["NNC1"])
+        cnnc2 = np.array(coag["NNC2"])
+        cnnct = np.array(coa["TRANNNC"])
+        coag = OpmGrid(dic["write"] + ".EGRID")
     rnnc1 = np.array(temp["NNC1"])
     rnnc2 = np.array(temp["NNC2"])
     rnnct = np.array(dic["ini"]["TRANNNC"])
-    cnnc1 = np.array(coag["NNC1"])
-    cnnc2 = np.array(coag["NNC2"])
-    cnnct = np.array(coa["TRANNNC"])
-    coag = OpmGrid(dic["write"] + ".EGRID")
-    refpv = np.array(dic["ini"]["PORV"])
     coapv = np.array(coa["PORV"])
-    coa_dz = np.zeros(len(coapv))
-    ref_dz = np.zeros(len(refpv))
     dic["coa_tranx"] = np.zeros(len(coapv))
     dic["coa_trany"] = np.zeros(len(coapv))
     dic["coa_tranx"][coapv > 0] = np.array(coa["TRANX"])
     dic["coa_trany"][coapv > 0] = np.array(coa["TRANY"])
-    coa_dz[coapv > 0] = np.array(coa["DZ"])
-    ref_dz[refpv > 0] = np.array(dic["ini"]["DZ"])
-    dic["coa_editnnc"] = []
-    dic["nnct"] = [[[] for _ in range(dic["yn"])] for _ in range(dic["xn"])]
-    indel = []
     for r1, r2, trn in zip(rnnc1, rnnc2, rnnct):
         rijk1 = dic["grid"].ijk_from_global_index(r1 - 1)
         rijk2 = dic["grid"].ijk_from_global_index(r2 - 1)
@@ -690,61 +691,61 @@ def handle_nnc_trans(dic):
                     + (dic["jc"][rijk1[1] + 1] - 1) * dic["nx"]
                     + (dic["kc"][rijk1[2] + 1] - 1) * dic["nx"] * dic["ny"]
                 )
-                if coa_dz[ind] > 0:
-                    dic["coa_tranx"][ind] += trn * ref_dz[r1 - 1] / coa_dz[ind]
+                dic["coa_tranx"][ind] += trn
             elif rijk1[1] + 1 == rijk2[1]:
                 ind = (
                     (dic["ic"][rijk1[0] + 1] - 1)
                     + (dic["jc"][rijk1[1] + 1] - 1) * dic["nx"]
                     + (dic["kc"][rijk1[2] + 1] - 1) * dic["nx"] * dic["ny"]
                 )
-                if coa_dz[ind] > 0:
-                    dic["coa_trany"][ind] += trn * ref_dz[r1 - 1] / coa_dz[ind]
+                dic["coa_trany"][ind] += trn
             elif rijk1[0] == rijk2[0] + 1:
                 ind = (
                     (dic["ic"][rijk2[0] + 1] - 1)
                     + (dic["jc"][rijk2[1] + 1] - 1) * dic["nx"]
                     + (dic["kc"][rijk2[2] + 1] - 1) * dic["nx"] * dic["ny"]
                 )
-                if coa_dz[ind] > 0:
-                    dic["coa_tranx"][ind] += trn * ref_dz[r2 - 1] / coa_dz[ind]
+                dic["coa_tranx"][ind] += trn
             elif rijk1[1] == rijk2[1] + 1:
                 ind = (
                     (dic["ic"][rijk2[0] + 1] - 1)
                     + (dic["jc"][rijk2[1] + 1] - 1) * dic["nx"]
                     + (dic["kc"][rijk2[2] + 1] - 1) * dic["nx"] * dic["ny"]
                 )
-                if coa_dz[ind] > 0:
-                    dic["coa_trany"][ind] += trn * ref_dz[r2 - 1] / coa_dz[ind]
-        else:
+                dic["coa_trany"][ind] += trn
+        elif editnnc:
             dic["nnct"][rijk1[0]][rijk1[1]].append([rijk1[2], rijk2[2], trn])
-    print("Processing the transmissibilities")
-    with alive_bar(len(cnnc1)) as bar_animation:
-        for n, (n1, n2) in enumerate(zip(cnnc1, cnnc2)):
-            bar_animation()
-            ijk1 = coag.ijk_from_global_index(n1 - 1)
-            ijk2 = coag.ijk_from_global_index(n2 - 1)
-            fip1, fip2 = ijk1[2] + 1, ijk2[2] + 1
-            rtran, found, indel = 0, 0, []
-            if fip1 != fip2 and (ijk1[0] != ijk2[0] or ijk1[1] != ijk2[1]):
-                for i, val in enumerate(dic["nnct"][ijk1[0]][ijk1[1]]):
-                    rfip1 = dic["kc"][val[0] + 1]
-                    rfip2 = dic["kc"][val[1] + 1]
-                    if fip1 == rfip1 and fip2 == rfip2:
-                        rtran += val[2]
-                        found = 1
-                        indel.append(i)
-                for ind in indel[::-1]:
-                    del dic["nnct"][ijk1[0]][ijk1[1]][ind]
-                if found == 1:
-                    mult = rtran / cnnct[n]
-                else:
-                    mult = 0
-                dic["coa_editnnc"].append(
-                    f"{ijk1[0]+1} {ijk1[1]+1} {ijk1[2]+1} {ijk2[0]+1} "
-                    f"{ijk2[1]+1} {ijk2[2]+1} {mult} /"
-                )
-    for name in ["tranx", "trany", "editnnc"]:
+    if editnnc:
+        print("Processing the transmissibilities")
+        with alive_bar(len(cnnc1)) as bar_animation:
+            for n, (n1, n2) in enumerate(zip(cnnc1, cnnc2)):
+                bar_animation()
+                ijk1 = coag.ijk_from_global_index(n1 - 1)
+                ijk2 = coag.ijk_from_global_index(n2 - 1)
+                fip1, fip2 = ijk1[2] + 1, ijk2[2] + 1
+                rtran, found, indel = 0, 0, []
+                if fip1 != fip2 and (ijk1[0] != ijk2[0] or ijk1[1] != ijk2[1]):
+                    for i, val in enumerate(dic["nnct"][ijk1[0]][ijk1[1]]):
+                        rfip1 = dic["kc"][val[0] + 1]
+                        rfip2 = dic["kc"][val[1] + 1]
+                        if fip1 == rfip1 and fip2 == rfip2:
+                            rtran += val[2]
+                            found = 1
+                            indel.append(i)
+                    for ind in indel[::-1]:
+                        del dic["nnct"][ijk1[0]][ijk1[1]][ind]
+                    if found == 1:
+                        mult = rtran / cnnct[n]
+                    else:
+                        mult = 0
+                    dic["coa_editnnc"].append(
+                        f"{ijk1[0]+1} {ijk1[1]+1} {ijk1[2]+1} {ijk2[0]+1} "
+                        f"{ijk2[1]+1} {ijk2[2]+1} {mult} /"
+                    )
+    cases = ["tranx", "trany"]
+    if editnnc:
+        cases += ["editnnc"]
+    for name in cases:
         if name == "editnnc":
             dic[f"coa_{name}"] = [f"{val}\n" for val in dic[f"coa_{name}"]]
         else:
@@ -853,7 +854,9 @@ def write_props(dic):
         for name in names:
             bar_animation()
             dic[f"{name}_c"] = compact_format(dic[f"{name}_c"])
-            if "*" in dic[f"{name}_c"][0]:
+            if "*" in dic[f"{name}_c"][0] and not (
+                dic["trans"] > 0 and name in ["tranx", "trany"]
+            ):
                 if int(dic[f"{name}_c"][0].split("*")[0]) == dic["ntot"]:
                     whr = dic["lol"].index(f"'{dic['label']}{name.upper()}.INC' /\n")
                     del dic["lol"][whr]

tests/test_7_transmissibilities.py

@@ -15,20 +15,12 @@
 
 
 def test_transmissibilities(flow):
-    """See examples/decks/MODEL4.DATA"""
+    """See examples/decks/MODEL4.DATA and MODEL5.DATA"""
     if not os.path.exists(f"{testpth}/output"):
         os.system(f"mkdir {testpth}/output")
-    subprocess.run(
-        [
-            flow,
-            f"{mainpth}/examples/decks/MODEL4.DATA",
-            f"--output-dir={testpth}/output/transmissibilities",
-        ],
-        check=True,
-    )
     values = [
-        [523.75178, 706.33996, 5008.1327, 106.21933],
-        [523.75178, 706.33996, 0, 106.21933],
+        [696.14886, 774.826, 5008.1323, 106.21955],
+        [696.14886, 774.826, 0, 106.21955],
     ]
     for i, sub in enumerate(["1", "2"]):
         subprocess.run(
@@ -58,5 +50,38 @@ def test_transmissibilities(flow):
         init = OpmFile(f"{testpth}/output/transmissibilities/COARSER{sub}.INIT")
         for j, n in enumerate(["X", "Y", "Z", "NNC"]):
             assert (
-                abs(np.sum(np.array(init[f"TRAN{n}"])) - values[i][j]) < 1e-2
+                abs(np.sum(np.array(init[f"TRAN{n}"])) - values[i][j]) < 1e-5
             ), f"Issue in TRAN{n} with -t {i}"
+        subprocess.run(
+            [
+                "pycopm",
+                "-f",
+                flow,
+                "-o",
+                f"{testpth}/output/transmissibilities",
+                "-i",
+                f"{mainpth}/examples/decks/MODEL5.DATA",
+                "-z",
+                "1:2",
+                "-w",
+                "COARSER_MODEL5",
+                "-l",
+                "CM5",
+                "-t",
+                "2",
+                "-a",
+                "max",
+                "-m",
+                "all",
+            ],
+            check=True,
+        )
+        ref = "MODEL5_PREP_PYCOPM_DRYRUN.INIT"
+        init_ref = OpmFile(f"{testpth}/output/transmissibilities/{ref}")
+        init = OpmFile(f"{testpth}/output/transmissibilities/COARSER_MODEL5.INIT")
+        assert (
+            abs(init["TRANX"][0] - 2 * init_ref["TRANX"][0]) < 1e-5
+        ), "Issue in TRANX with MODEL5 (NNCTRAN[0])"
+        assert (
+            abs(init["TRANX"][-2] - 2 * init_ref["TRANX"][0]) < 1e-5
+        ), "Issue in TRANX with MODEL5 (NNCTRAN[-2])"