Merge pull request #22 from daavid00/reportData

Dense data from all simulation grid cells

examples/finner_grids/spe11a.txt

@@ -17,10 +17,10 @@ cartesian         #Type of grid (cartesian, tensor, or corner-point)
 0 0               #Elevation of the parabola [m] (for spe11c)
 
 """Set the saturation functions"""
-(max(0, (s_w - swi) / (1 - swi))) ** 1.5                                                        #Wetting rel perm saturation function [-]
-(max(0, (1 - s_w - sni) / (1 - sni))) ** 1.5                                                    #Non-wetting rel perm saturation function [-]
-penmax * math.erf(pen * ((s_w-swi) / (1.-swi)) ** (-(1.0 / 1.5)) * math.pi**0.5 / (penmax * 2)) #Capillary pressure saturation function [Pa]
-(np.exp(np.flip(np.linspace(0, 5.0, npoints))) - 1) / (np.exp(5.0) - 1)                         #Points to evaluate the saturation functions (s_w) [-]
+(max(0, (s_w - swi) / (1 - swi))) ** 2                                                        #Wetting rel perm saturation function [-]
+(max(0, (1 - s_w - sni) / (1 - sni))) ** 2                                                    #Non-wetting rel perm saturation function [-]
+penmax * math.erf(pen * ((s_w-swi) / (1.-swi)) ** (-(1.0 / 2)) * math.pi**0.5 / (penmax * 2)) #Capillary pressure saturation function [Pa]
+(np.exp(np.flip(np.linspace(0, 5.0, npoints))) - 1) / (np.exp(5.0) - 1)                       #Points to evaluate the saturation functions (s_w) [-]
 
 """Properties sat functions"""
 """swi [-], sni [-], pen [Pa], penmax [Pa], npoints [-]"""

examples/hello_world/spe11a.txt

@@ -17,20 +17,20 @@ tensor            #Type of grid (cartesian, tensor, or corner-point)
 0 0               #Elevation of the parabola [m] (for spe11c) 
 
 """Set the saturation functions"""
-(max(0, (s_w - swi) / (1 - swi))) ** 1.5                                                        #Wetting rel perm saturation function [-]
-(max(0, (1 - s_w - sni) / (1 - sni))) ** 1.5                                                    #Non-wetting rel perm saturation function [-]
-penmax * math.erf(pen * ((s_w-swi) / (1.-swi)) ** (-(1.0 / 1.5)) * math.pi**0.5 / (penmax * 2)) #Capillary pressure saturation function [Pa]
-(np.exp(np.flip(np.linspace(0, 5.0, npoints))) - 1) / (np.exp(5.0) - 1)                         #Points to evaluate the saturation functions (s_w) [-]
+(max(0, (s_w - swi) / (1 - swi))) ** 2                                                        #Wetting rel perm saturation function [-]
+(max(0, (1 - s_w - sni) / (1 - sni))) ** 2                                                    #Non-wetting rel perm saturation function [-]
+penmax * math.erf(pen * ((s_w-swi) / (1.-swi)) ** (-(1.0 / 2)) * math.pi**0.5 / (penmax * 2)) #Capillary pressure saturation function [Pa]
+(np.exp(np.flip(np.linspace(0, 5.0, npoints))) - 1) / (np.exp(5.0) - 1)                       #Points to evaluate the saturation functions (s_w) [-]
 
 """Properties sat functions"""
 """swi [-], sni [-], pen [Pa], penmax [Pa], npoints [-]"""
-SWI1 0.32 SNI1 0.1 PEN1 1500 PENMAX1 9.5e4 NPOINTS1 100
-SWI2 0.14 SNI2 0.1 PEN2  300 PENMAX2 9.5e4 NPOINTS2 100
-SWI3 0.12 SNI3 0.1 PEN3  100 PENMAX3 9.5e4 NPOINTS3 100
-SWI4 0.12 SNI4 0.1 PEN4   25 PENMAX4 9.5e4 NPOINTS4 100
-SWI5 0.12 SNI5 0.1 PEN5   10 PENMAX5 9.5e4 NPOINTS5 100
-SWI6 0.10 SNI6 0.1 PEN6    1 PENMAX6 9.5e4 NPOINTS6 100
-SWI7    0 SNI7   0 PEN7    0 PENMAX7 9.5e4 NPOINTS7   2
+SWI1 0.32 SNI1 0.1 PEN1 1500 PENMAX1 2500 NPOINTS1 100
+SWI2 0.14 SNI2 0.1 PEN2  300 PENMAX2 2500 NPOINTS2 100
+SWI3 0.12 SNI3 0.1 PEN3  100 PENMAX3 2500 NPOINTS3 100
+SWI4 0.12 SNI4 0.1 PEN4   25 PENMAX4 2500 NPOINTS4 100
+SWI5 0.12 SNI5 0.1 PEN5   10 PENMAX5 2500 NPOINTS5 100
+SWI6 0.10 SNI6 0.1 PEN6    1 PENMAX6 2500 NPOINTS6 100
+SWI7    0 SNI7   0 PEN7    0 PENMAX7 2500 NPOINTS7   2
 
 """Properties rock"""
 """K [mD], phi [-]"""

src/pyopmspe11/utils/inputvalues.py

@@ -28,13 +28,13 @@ def process_input(dic, in_file):
             lol.append(row)
     dic = readthefirstpart(lol, dic)
     if dic["spe11"] == "spe11a":
-        dic["sensors"] = [[1.5, 0.0, 0.5], [1.7, 0.0, 1.1]]
-        dic["boxa"] = [[1.1, 0.0, 0.0], [2.8, 1.0, 0.6]]
-        dic["boxb"] = [[0.0, 0.0, 0.6], [1.1, 1.0, 1.2]]
-        dic["boxc"] = [[1.1, 0.0, 0.1], [2.6, 1.0, 0.4]]
+        dic["sensors"] = [[1.5, 0.005, 0.5], [1.7, 0.005, 1.1]]
+        dic["boxa"] = [[1.1, 0.0, 0.0], [2.8, 0.01, 0.6]]
+        dic["boxb"] = [[0.0, 0.0, 0.6], [1.1, 0.01, 1.2]]
+        dic["boxc"] = [[1.1, 0.0, 0.1], [2.6, 0.01, 0.4]]
         dic["time"] = 3600.0  # hour to seconds
     elif dic["spe11"] == "spe11b":
-        dic["sensors"] = [[4500.0, 0.0, 500], [5100.0, 0.0, 1100.0]]
+        dic["sensors"] = [[4500.0, 0.5, 500], [5100.0, 0.5, 1100.0]]
         dic["boxa"] = [[3300.0, 0.0, 0.0], [8300.0, 1.0, 600.0]]
         dic["boxb"] = [[100.0, 0.0, 600.0], [3300.0, 1.0, 1200.0]]
         dic["boxc"] = [[3300.0, 0.0, 100.0], [7800.0, 1.0, 400.0]]

src/pyopmspe11/visualization/data.py

@@ -736,25 +736,35 @@ def dense_data(dic):
     for i, j, k in zip(["x", "y", "z"], dic["dims"], dic["nxyz"]):
         temp = np.linspace(0, j, k + 1)
         dic[f"ref{i}cent"] = 0.5 * (temp[1:] + temp[:-1])
-    ind, dic["cell_ind"] = 0, [0] * dic["nocellsr"]
-    # This is very slow, it needs to be fixed
-    for k in np.flip(dic["refzcent"]):
+        dic[f"ref{i}grid"] = []
+    for k in dic["refzcent"]:
         for j in dic["refycent"]:
             for i in dic["refxcent"]:
-                dic["cell_ind"][ind] = pd.Series(
-                    (dic["simxcent"] - i) ** 2
-                    + (dic["simycent"] - j) ** 2
-                    + (dic["simzcent"] - k) ** 2
-                ).argmin()
-                ind += 1
+                dic["refxgrid"].append(i)
+                dic["refygrid"].append(j)
+                dic["refzgrid"].append(k)
+    for i in ["x", "y", "z"]:
+        dic[f"ref{i}grid"] = np.array(dic[f"ref{i}grid"])
+    ind, dic["cell_ind"] = 0, [0] * dic["nocellst"]
+    # This is very slow, it needs to be fixed
+    for i, j, k in zip(dic["simxcent"], dic["simycent"], np.flip(dic["simzcent"])):
+        dic["cell_ind"][ind] = pd.Series(
+            (dic["refxgrid"] - i) ** 2
+            + (dic["refygrid"] - j) ** 2
+            + (dic["refzgrid"] - k) ** 2
+        ).argmin()
+        ind += 1
+    dic["porv_array"] = np.array(dic["porv"])
+    if max(dic["satnum"]) < 7:
+        dic = handle_inactive_mapping(dic)
     names = ["pressure", "sgas", "xco2", "xh20", "gden", "wden", "tco2"]
     if dic["case"] != "spe11a":
         names += ["temp"]
     for i in range(n_t + 1):
         print(f"Processing dense data {i+1} out of {n_t + 1}")
         t_n = i * d_s + dic["t1_rst"]
         for name in names:
-            dic[f"{name}_array"] = np.array([np.nan] * dic["nocellst"])
+            dic[f"{name}_array"] = np.array([0.0] * dic["nocellst"])
         co2_g = [
             sga * rho * por
             for (sga, rho, por) in zip(
@@ -807,12 +817,48 @@ def dense_data(dic):
                 h2o / (h2o + co2) if (h2o + co2) > 0 else 0.0
                 for (h2o, co2) in zip(h2o_v, co2_g)
             ]
-        for name in names:
-            dic[f"{name}_refg"] = [np.nan] * dic["nocellsr"]
-            dic[f"{name}_refg"] = dic[f"{name}_array"][dic["cell_ind"]]
+        dic = map_to_report_grid(dic, names)
         write_dense_data(dic, i)
 
 
+def handle_inactive_mapping(dic):
+    """Set to inf the inactive grid centers in the reporting grid"""
+    var_array = np.array([np.nan] * dic["nocellst"])
+    var_array[dic["actind"]] = dic["sgas"][0]
+    for ii in range(dic["nocellsr"]):
+        inds = [iii == ii for iii in dic["cell_ind"]]
+        if np.isnan(sum(var_array[inds])):
+            dic["refxgrid"][ii] = np.inf
+            dic["refygrid"][ii] = np.inf
+            dic["refzgrid"][ii] = np.inf
+    ind, dic["cell_ind"] = 0, [0] * dic["nocellst"]
+    for i, j, k in zip(dic["simxcent"], dic["simycent"], np.flip(dic["simzcent"])):
+        dic["cell_ind"][ind] = pd.Series(
+            (dic["refxgrid"] - i) ** 2
+            + (dic["refygrid"] - j) ** 2
+            + (dic["refzgrid"] - k) ** 2
+        ).argmin()
+        ind += 1
+    return dic
+
+
+def map_to_report_grid(dic, names):
+    """Map the simulation grid to the reporting grid"""
+    for name in names:
+        dic[f"{name}_refg"] = [np.nan] * dic["nocellsr"]
+        for ii in range(dic["nocellsr"]):
+            inds = [iii == ii for iii in dic["cell_ind"]]
+            p_v = sum(dic["porv_array"][inds])
+            if p_v > 0:
+                if name == "tco2":
+                    dic[f"{name}_refg"][ii] = sum(dic[f"{name}_array"][inds])
+                else:
+                    dic[f"{name}_refg"][ii] = (
+                        sum(dic[f"{name}_array"][inds] * dic["porv_array"][inds]) / p_v
+                    )
+    return dic
+
+
 def write_dense_data(dic, i):
     """Map the quantities to the cells"""
     if dic["case"] == "spe11a":

tests/configs/input.txt

@@ -40,7 +40,7 @@ PERM3 202.650 PORO3 0.20 THCONR3 1.25
 PERM4 506.625 PORO4 0.20 THCONR4 1.25
 PERM5 1013.25 PORO5 0.25 THCONR5 0.92
 PERM6 2026.50 PORO6 0.35 THCONR6 0.26
-PERM7       0 PORO7    0 THCONR7 2.00
+PERM7    1e-5 PORO7 1e-6 THCONR7 2.00
 
 """Wells radius and position"""
 """radius (0 to use the SOURCE keyword instead of well keywords, this requires a Flow version newer than 23-01-2024), x, y, and z position [m] (final positions as well for spe11c)"""

tests/configs/spe11a.txt

@@ -17,10 +17,10 @@ cartesian         #Type of grid (cartesian or corner-point)
 0 0               #Elevation of the parabola [m] (for spe11c) 
 
 """Set the saturation functions"""
-(max(0, (s_w - swi) / (1 - swi))) ** 1.5                                                        #Wetting rel perm saturation function [-]
-(max(0, (1 - s_w - sni) / (1 - sni))) ** 1.5                                                    #Non-wetting rel perm saturation function [-]
-penmax * math.erf(pen * ((s_w-swi) / (1.-swi)) ** (-(1.0 / 1.5)) * math.pi**0.5 / (penmax * 2)) #Capillary pressure saturation function [Pa]
-(np.exp(np.flip(np.linspace(0, 5.0, npoints))) - 1) / (np.exp(5.0) - 1)                         #Points to evaluate the saturation functions (s_w) [-]
+(max(0, (s_w - swi) / (1 - swi))) ** 2                                                        #Wetting rel perm saturation function [-]
+(max(0, (1 - s_w - sni) / (1 - sni))) ** 2                                                    #Non-wetting rel perm saturation function [-]
+penmax * math.erf(pen * ((s_w-swi) / (1.-swi)) ** (-(1.0 / 2)) * math.pi**0.5 / (penmax * 2)) #Capillary pressure saturation function [Pa]
+(np.exp(np.flip(np.linspace(0, 5.0, npoints))) - 1) / (np.exp(5.0) - 1)                       #Points to evaluate the saturation functions (s_w) [-]
 
 """Properties sat functions"""
 """swi [-], sni [-], pen [Pa], penmax [Pa], npoints [-]"""

tests/configs/spe11c.txt

@@ -40,7 +40,7 @@ PERM3 202.650 PORO3 0.20 THCONR3 1.25
 PERM4 506.625 PORO4 0.20 THCONR4 1.25
 PERM5 1013.25 PORO5 0.25 THCONR5 0.92
 PERM6 2026.50 PORO6 0.35 THCONR6 0.26
-PERM7       0 PORO7    0 THCONR7 2.00
+PERM7    1e-5 PORO7 1e-6 THCONR7 2.00
 
 """Wells radius and position"""
 """radius (0 to use the SOURCE keyword instead of well keywords, this requires a Flow version newer than 23-01-2024), x, y, and z position [m] (final positions as well for spe11c)"""
@@ -49,5 +49,4 @@ PERM7       0 PORO7    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), ..."""
- 5  5 1 1  0 10 1  0 10
-25  5 1 1 50 10 1 50 10
+25 5 1 1 50 10 1 50 10

tests/test_2_spe11c.py

@@ -23,7 +23,7 @@ def test_spe11b():
             "-g",
             "all",
             "-r",
-            "50,5,25",
+            "24,1,12",
             "-t",
             "5",
         ],