Merge pull request #45 from totto82/fixPVB

Include porosity in the pore volume boundary layer

src/pyopmspe11/utils/mapproperties.py

@@ -173,7 +173,9 @@ def structured_handling_spe11bc(dic):
             dic["satnum"].append(dic["ids_gmsh"][fgl][0])
             dic = boxes(dic, dic["xmx_center"][i], z_c, i, dic["satnum"][-1])
             dic["permx"].append(dic["rock"][int(dic["ids_gmsh"][fgl][0]) - 1][0])
-            dic["poro"].append(dic["rock"][int(dic["ids_gmsh"][fgl][0]) - 1][1])
+            poro = dic["rock"][int(dic["ids_gmsh"][fgl][0]) - 1][1]
+            dic["poro"].append(poro)
+            pv = float(poro) * (dic["pvAdded"] + dic["widthBuffer"])
             dic["thconr"].append(
                 f"{dic['rockCond'][int(dic['ids_gmsh'][fgl][0])-1][0]}"
             )
@@ -184,13 +186,13 @@ def structured_handling_spe11bc(dic):
                 int(dic["ids_gmsh"][fgl][0]) != 1 and int(dic["ids_gmsh"][fgl][0]) != 7
             ):
                 dic["porv"].append(
-                    f"PORV {dic['pvAdded']*dic['dy'][0]*dic['dz'][k]} 1 1 1 1 {k+1} {k+1} /"
+                    f"PORV {pv*dic['dy'][0]*dic['dz'][k]} 1 1 1 1 {k+1} {k+1} /"
                 )
             elif i == dic["noCells"][0] - 1 and (
                 int(dic["ids_gmsh"][fgl][0]) != 1 and int(dic["ids_gmsh"][fgl][0]) != 7
             ):
                 dic["porv"].append(
-                    f"PORV {dic['pvAdded']*dic['dy'][0]*dic['dz'][k]} {dic['noCells'][0]} "
+                    f"PORV {pv*dic['dy'][0]*dic['dz'][k]} {dic['noCells'][0]} "
                     + f"{dic['noCells'][0]} 1 1 {k+1} {k+1} /"
                 )
             centers.append(
@@ -235,15 +237,15 @@ def structured_handling_spe11bc(dic):
                     and int(dic["satnum"][-dic["noCells"][0] + i_i]) != 7
                 ):
                     dic["porv"].append(
-                        f"PORV {dic['pvAdded']*dic['dy'][j+1]*dic['dz'][k]} 1 1 "
+                        f"PORV {pv*dic['dy'][j+1]*dic['dz'][k]} 1 1 "
                         + f"{j+2} {j+2} {k+1} {k+1} /"
                     )
                 elif i_i == dic["noCells"][0] - 1 and (
                     int(dic["satnum"][-dic["noCells"][0] + i_i]) != 1
                     and int(dic["satnum"][-dic["noCells"][0] + i_i]) != 7
                 ):
                     dic["porv"].append(
-                        f"PORV {dic['pvAdded']*dic['dy'][j+1]*dic['dz'][k]} {dic['noCells'][0]} "
+                        f"PORV {pv*dic['dy'][j+1]*dic['dz'][k]} {dic['noCells'][0]} "
                         + f"{dic['noCells'][0]} {j+2} {j+2} {k+1} {k+1} /"
                     )
     dic["pop1"] = pd.Series(sensor1).argmin()
@@ -435,21 +437,23 @@ def corner_point_handling_spe11bc(dic):
         dic["satnum"].append(dic["ids_gmsh"][fgl][0])
         dic = boxes(dic, dic["xyz"][0], z_c, dic["ijk"][0], dic["satnum"][-1])
         dic["permx"].append(dic["rock"][int(dic["ids_gmsh"][fgl][0]) - 1][0])
-        dic["poro"].append(dic["rock"][int(dic["ids_gmsh"][fgl][0]) - 1][1])
+        poro = dic["rock"][int(dic["ids_gmsh"][fgl][0]) - 1][1]
+        dic["poro"].append(poro)
+        pv = float(poro) * (dic["pvAdded"] + dic["widthBuffer"])
         dic["disperc"].append(f"{dic['dispersion'][int(dic['ids_gmsh'][fgl][0])-1]}")
         dic["thconr"].append(f"{dic['rockCond'][int(dic['ids_gmsh'][fgl][0])-1][0]}")
         if dic["ijk"][0] == 0 and (
             int(dic["ids_gmsh"][fgl][0]) != 1 and int(dic["ids_gmsh"][fgl][0]) != 7
         ):
             dic["porv"].append(
-                f"PORV {dic['pvAdded']*dic['d_y'][0]*dic['d_z'][i]} 1 1 1 1 "
+                f"PORV { pv*dic['d_y'][0]*dic['d_z'][i]} 1 1 1 1 "
                 + f"{dic['ijk'][2]+1} {dic['ijk'][2]+1} /"
             )
         elif dic["ijk"][0] == dic["noCells"][0] - 1 and (
             int(dic["ids_gmsh"][fgl][0]) != 1 and int(dic["ids_gmsh"][fgl][0]) != 7
         ):
             dic["porv"].append(
-                f"PORV {dic['pvAdded']*dic['d_y'][0]*dic['d_z'][i]} {dic['noCells'][0]} "
+                f"PORV {pv*dic['d_y'][0]*dic['d_z'][i]} {dic['noCells'][0]} "
                 + f"{dic['noCells'][0]} 1 1 {dic['ijk'][2]+1} {dic['ijk'][2]+1} /"
             )
         centers.append(str([dic["xyz"][0], dic["ymy_center"][0], dic["xyz"][2]])[1:-1])
@@ -476,15 +480,15 @@ def corner_point_handling_spe11bc(dic):
                         dic["d_zl"] = dic["d_z"][-dic["noCells"][0] + 1 + i]
                         dic["porv"].append(
                             "PORV "
-                            + f"{dic['pvAdded']*dic['d_y'][j+1]*dic['d_zl']} 1 1 "
+                            + f"{pv*dic['d_y'][j+1]*dic['d_zl']} 1 1 "
                             + f"{j+2} {j+2} {dic['ijk'][2]+1} {dic['ijk'][2]+1} /"
                         )
                     elif i_i == dic["noCells"][0] - 1 and (
                         int(dic["satnum"][-dic["noCells"][0] + i_i]) != 1
                         and int(dic["satnum"][-dic["noCells"][0] + i_i]) != 7
                     ):
                         dic["porv"].append(
-                            f"PORV {dic['pvAdded']*dic['d_y'][j+1]*dic['d_z'][i]} "
+                            f"PORV {pv*dic['d_y'][j+1]*dic['d_z'][i]} "
                             + f"{dic['noCells'][0]} {dic['noCells'][0]} {j+2} {j+2} "
                             + f"{dic['ijk'][2]+1} {dic['ijk'][2]+1} /"
                         )

tests/configs/input.txt

@@ -49,4 +49,4 @@ PERM7    1e-5 PORO7 1e-6 DISP7 0 THCONR7 2.00
 
 """Define the injection values ([hours] for spe11a; [years] for spe11b/c)""" 
 """injection time, time step size to write results, maximum solver time step, injected fluid (0 water, 1 co2) (well1), injection rate [kg/s] (well1), temperature [C] (well1), injected fluid (0 water, 1 co2) (well2), ..."""
-25 5 1 1 0.035 10 1 0 10
+25 5 .1 1 0.035 10 1 0 10