Merge pull request #34 from daavid00/dataReport

Improving the performance to generate the dense data

docs/_sources/examples.rst.txt

@@ -64,14 +64,29 @@ with grid size of the same order as requested for reporting the spatial data in
 
 .. code-block:: bash
 
-    pyopmspe11 -i spe11b.txt -o spe11b_finner -m all -g all -r 840,1,120 -t 5
+    pyopmspe11 -i spe11b.txt -o spe11b_current -m all -g all -r 840,1,120 -t 5 -w 0.1
 
 the following are some of the generated figures:
 
 .. figure:: figs/spe11b_sparse_data_finner.png
+
+    Sparse data.
+
 .. figure:: figs/spe11b_performance_finner.png
-.. figure:: figs/spe11b_xco2_2Dmaps_finner.png
 
-.. tip::
-    By executing flow --help you get an overview of the available flags in the flow simulator to improve/fix convergence issues 
-    (i.e., by setting the flag --linear-solver=cprw to change the linear solver).
+    Performance data.
+
+.. figure:: figs/spe11b_performance_detailed_finner.png
+
+    Detailed performance data.
+
+.. figure:: figs/spe11b_wden_2Dmaps_finner.png
+
+    Dense data for the water density.
+
+.. warning::
+    The inputs in this configuration file were set to show the current functionality of pyopmspe11 to generate the benchmark data, i.e.,
+    the tolerances are not tight and no tunning was performed. In consequence, the total mass in the simulation (ca. 7.75e7 kg) is less
+    than what is specified in the benchmark (ca. 8.27e7 kg). By executing flow --help you get an overview of the available flags in the
+    flow simulator to improve/fix convergence and mass issues (i.e., by setting the flag --linear-solver=cprw to change the linear solver,
+    by tightening the mb tolerances (--tolerance-mb), etc.).

docs/examples.html

@@ -135,23 +135,41 @@ <h2>Presentation<a class="headerlink" href="#presentation" title="Link to this h
 </div>
 <p>The <a class="reference external" href="https://github.com/OPM/pyopmspe11/blob/main/examples">examples/finner_grids</a> folder contains configuration files
 with grid size of the same order as requested for reporting the spatial data in the benchmark, as well as the required injection schedules. For example, by executing:</p>
-<div class="highlight-bash notranslate"><div class="highlight"><pre><span></span>pyopmspe11<span class="w"> </span>-i<span class="w"> </span>spe11b.txt<span class="w"> </span>-o<span class="w"> </span>spe11b_finner<span class="w"> </span>-m<span class="w"> </span>all<span class="w"> </span>-g<span class="w"> </span>all<span class="w"> </span>-r<span class="w"> </span><span class="m">840</span>,1,120<span class="w"> </span>-t<span class="w"> </span><span class="m">5</span>
+<div class="highlight-bash notranslate"><div class="highlight"><pre><span></span>pyopmspe11<span class="w"> </span>-i<span class="w"> </span>spe11b.txt<span class="w"> </span>-o<span class="w"> </span>spe11b_current<span class="w"> </span>-m<span class="w"> </span>all<span class="w"> </span>-g<span class="w"> </span>all<span class="w"> </span>-r<span class="w"> </span><span class="m">840</span>,1,120<span class="w"> </span>-t<span class="w"> </span><span class="m">5</span><span class="w"> </span>-w<span class="w"> </span><span class="m">0</span>.1
 </pre></div>
 </div>
 <p>the following are some of the generated figures:</p>
-<figure class="align-default">
+<figure class="align-default" id="id1">
 <img alt="_images/spe11b_sparse_data_finner.png" src="_images/spe11b_sparse_data_finner.png" />
+<figcaption>
+<p><span class="caption-text">Sparse data.</span><a class="headerlink" href="#id1" title="Link to this image"></a></p>
+</figcaption>
 </figure>
-<figure class="align-default">
+<figure class="align-default" id="id2">
 <img alt="_images/spe11b_performance_finner.png" src="_images/spe11b_performance_finner.png" />
+<figcaption>
+<p><span class="caption-text">Performance data.</span><a class="headerlink" href="#id2" title="Link to this image"></a></p>
+</figcaption>
 </figure>
-<figure class="align-default">
-<img alt="_images/spe11b_xco2_2Dmaps_finner.png" src="_images/spe11b_xco2_2Dmaps_finner.png" />
+<figure class="align-default" id="id3">
+<img alt="_images/spe11b_performance_detailed_finner.png" src="_images/spe11b_performance_detailed_finner.png" />
+<figcaption>
+<p><span class="caption-text">Detailed performance data.</span><a class="headerlink" href="#id3" title="Link to this image"></a></p>
+</figcaption>
+</figure>
+<figure class="align-default" id="id4">
+<img alt="_images/spe11b_wden_2Dmaps_finner.png" src="_images/spe11b_wden_2Dmaps_finner.png" />
+<figcaption>
+<p><span class="caption-text">Dense data for the water density.</span><a class="headerlink" href="#id4" title="Link to this image"></a></p>
+</figcaption>
 </figure>
-<div class="admonition tip">
-<p class="admonition-title">Tip</p>
-<p>By executing flow –help you get an overview of the available flags in the flow simulator to improve/fix convergence issues
-(i.e., by setting the flag –linear-solver=cprw to change the linear solver).</p>
+<div class="admonition warning">
+<p class="admonition-title">Warning</p>
+<p>The inputs in this configuration file were set to show the current functionality of pyopmspe11 to generate the benchmark data, i.e.,
+the tolerances are not tight and no tunning was performed. In consequence, the total mass in the simulation (ca. 7.75e7 kg) is less
+than what is specified in the benchmark (ca. 8.27e7 kg). By executing flow –help you get an overview of the available flags in the
+flow simulator to improve/fix convergence and mass issues (i.e., by setting the flag –linear-solver=cprw to change the linear solver,
+by tightening the mb tolerances (–tolerance-mb), etc.).</p>
 </div>
 </section>
 </section>

docs/text/examples.rst

@@ -64,14 +64,29 @@ with grid size of the same order as requested for reporting the spatial data in
 
 .. code-block:: bash
 
-    pyopmspe11 -i spe11b.txt -o spe11b_finner -m all -g all -r 840,1,120 -t 5
+    pyopmspe11 -i spe11b.txt -o spe11b_current -m all -g all -r 840,1,120 -t 5 -w 0.1
 
 the following are some of the generated figures:
 
 .. figure:: figs/spe11b_sparse_data_finner.png
+
+    Sparse data.
+
 .. figure:: figs/spe11b_performance_finner.png
-.. figure:: figs/spe11b_xco2_2Dmaps_finner.png
 
-.. tip::
-    By executing flow --help you get an overview of the available flags in the flow simulator to improve/fix convergence issues 
-    (i.e., by setting the flag --linear-solver=cprw to change the linear solver).
+    Performance data.
+
+.. figure:: figs/spe11b_performance_detailed_finner.png
+
+    Detailed performance data.
+
+.. figure:: figs/spe11b_wden_2Dmaps_finner.png
+
+    Dense data for the water density.
+
+.. warning::
+    The inputs in this configuration file were set to show the current functionality of pyopmspe11 to generate the benchmark data, i.e.,
+    the tolerances are not tight and no tunning was performed. In consequence, the total mass in the simulation (ca. 7.75e7 kg) is less
+    than what is specified in the benchmark (ca. 8.27e7 kg). By executing flow --help you get an overview of the available flags in the
+    flow simulator to improve/fix convergence and mass issues (i.e., by setting the flag --linear-solver=cprw to change the linear solver,
+    by tightening the mb tolerances (--tolerance-mb), etc.).

examples/finner_grids/spe11b.txt

@@ -1,5 +1,5 @@
 """Set the full path to the flow executable and flags"""
-mpirun -np 70 flow --tolerance-mb=1e-7 --linear-solver=cprw --newton-min-iterations=1 --enable-tuning=true --enable-opm-rst-file=true --output-extra-convergence-info=steps,iterations
+mpirun -np 70 flow --linear-solver=cprw --newton-min-iterations=1 --enable-tuning=true --enable-opm-rst-file=true --output-extra-convergence-info=steps,iterations
 
 """Set the model parameters"""
 spe11b master     #Name of the spe case (spe11a, spe11b, or spe11c) and OPM Flow version (master or release)

src/pyopmspe11/utils/mapproperties.py

@@ -180,7 +180,11 @@ def structured_handling_spe11bc(dic):
                     f"PORV {dic['pvAdded']*dic['dy'][0]*dic['dz'][k]} {dic['noCells'][0]} "
                     + f"{dic['noCells'][0]} 1 1 {k+1} {k+1} /"
                 )
-            centers.append(str([dic["xmx_center"][i], dic["ymy_center"][0], z_c])[1:-1])
+            centers.append(
+                str([dic["xmx_center"][i], dic["ymy_center"][0], dic["zmz_center"][k]])[
+                    1:-1
+                ]
+            )
         for j in range(dic["noCells"][1] - 1):
             for names in ["satnum", "poro", "permx", "disperc", "thconr"]:
                 dic[f"{names}"].extend(dic[f"{names}"][-dic["noCells"][0] :])
@@ -207,9 +211,6 @@ def structured_handling_spe11bc(dic):
                     i_i,
                     dic["satnum"][-dic["noCells"][0] + i_i],
                 )
-                centers.append(
-                    str([dic["xmx_center"][i_i], dic["ymy_center"][j + 1], z_c])[1:-1]
-                )
                 if i_i == 0 and (
                     int(dic["satnum"][-dic["noCells"][0] + i_i]) != 1
                     and int(dic["satnum"][-dic["noCells"][0] + i_i]) != 7
@@ -399,7 +400,7 @@ def corner_point_handling_spe11bc(dic):
                 f"PORV {dic['pvAdded']*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], z_c])[1:-1])
+        centers.append(str([dic["xyz"][0], dic["ymy_center"][0], dic["xyz"][2]])[1:-1])
         if dic["ijk"][0] > 0 and dic["ijk"][0] == dic["noCells"][0] - 1:
             for j in range(dic["noCells"][1] - 1):
                 for names in ["satnum", "poro", "permx", "disperc", "thconr"]:
@@ -415,9 +416,6 @@ def corner_point_handling_spe11bc(dic):
                         i_i,
                         dic["satnum"][-dic["noCells"][0] + i_i],
                     )
-                    centers.append(
-                        str([xtemp[i_i], dic["ymy_center"][j + 1], z_c])[1:-1]
-                    )
                     if i_i == 0 and (
                         int(dic["satnum"][-dic["noCells"][0] + i_i]) != 1
                         and int(dic["satnum"][-dic["noCells"][0] + i_i]) != 7
@@ -598,6 +596,9 @@ def positions(dic):
             dic = structured_handling_spe11a(dic)
         else:
             dic = structured_handling_spe11bc(dic)
+    np.savetxt(
+        f"{dic['exe']}/{dic['fol']}/deck/ycenters.txt", dic["ymy_center"], fmt="%.8E"
+    )
     return dic
 
 

src/pyopmspe11/visualization/data.py

@@ -181,6 +181,7 @@ def read_resdata(dig):
         dig["egrid"].ny,
         dig["egrid"].nz,
     ]
+    dig["noxz"] = dig["egrid"].nx * dig["egrid"].nz
     return dig
 
 
@@ -204,6 +205,7 @@ def read_opm(dig):
         dig["egrid"].dimension[1],
         dig["egrid"].dimension[2],
     ]
+    dig["noxz"] = dig["egrid"].dimension[0] * dig["egrid"].dimension[2]
     return dig
 
 
@@ -651,12 +653,11 @@ def max_xcw(dig, dil):
 
 def get_centers(dig, dil):
     """Centers from the simulation grid"""
-    for i in ["x", "y", "z"]:
-        dil[f"sim{i}cent"] = [0.0] * dig["nocellst"]
+    for i in ["x", "z"]:
+        dil[f"sim{i}cent"] = [0.0] * dig["noxz"]
     with open(f"{dig['path']}/deck/centers.txt", "r", encoding="utf8") as file:
         for j, row in enumerate(csv.reader(file)):
             dil["simxcent"][j] = float(row[0])
-            dil["simycent"][j] = float(row[1])
             dil["simzcent"][j] = dig["dims"][2] - float(row[2])
     if dig["use"] == "opm":
         dil["satnum"] = list(dig["init"]["SATNUM"])
@@ -666,7 +667,12 @@ def get_centers(dig, dil):
 
 
 def dense_data(dig):
-    """Write the quantities with the benchmark format"""
+    """
+    Write the quantities with the benchmark format for the dense data.
+    Still plenty of room to improve here the performance and memory
+    usage. However, this current implementation can handle cases with
+    ca. 15 000 000 cells in a reasonable time (ca. hours).
+    """
     dil = {"rstno": []}
     for time in dig["dense_t"]:
         dil["rstno"].append(dig["times"].index(time))
@@ -675,25 +681,24 @@ def dense_data(dig):
     for i, j, k in zip(["x", "y", "z"], dig["dims"], dig["nxyz"]):
         ind = np.linspace(0, j, k + 1)
         dil[f"ref{i}cent"] = 0.5 * (ind[1:] + ind[:-1])
-        dil[f"ref{i}grid"] = np.zeros(dig["nocellsr"])
+    dil["refxgrid"] = np.zeros(dig["nxyz"][0] * dig["nxyz"][2])
+    dil["refzgrid"] = np.zeros(dig["nxyz"][0] * dig["nxyz"][2])
     ind = 0
     for k in dil["refzcent"]:
-        for j in dil["refycent"]:
-            for i in dil["refxcent"]:
-                dil["refxgrid"][ind] = i
-                dil["refygrid"][ind] = j
-                dil["refzgrid"][ind] = k
-                ind += 1
-    ind, dil["cell_ind"] = 0, np.zeros(dig["nocellst"], dtype=int)
-    for i, j, k in zip(dil["simxcent"], dil["simycent"], dil["simzcent"]):
+        for i in dil["refxcent"]:
+            dil["refxgrid"][ind] = i
+            dil["refzgrid"][ind] = k
+            ind += 1
+    ind, dil["cell_ind"] = 0, np.zeros(dig["noxz"], dtype=int)
+    for i, k in zip(dil["simxcent"], dil["simzcent"]):
         dil["cell_ind"][ind] = pd.Series(
-            (dil["refxgrid"] - i) ** 2
-            + (dil["refygrid"] - j) ** 2
-            + (dil["refzgrid"] - k) ** 2
+            np.abs(dil["refxgrid"] - i) + np.abs(dil["refzgrid"] - k)
         ).argmin()
         ind += 1
-    if max(dil["satnum"]) < 7:
+    if max(dil["satnum"]) < 7 and dig["case"] == "spe11a":
         dil = handle_inactive_mapping(dig, dil)
+    if dig["case"] == "spe11c":
+        dil = handle_yaxis_mapping(dig, dil)
     if dig["mode"] == "all" or dig["mode"][:5] == "dense":
         names = ["pressure", "sgas", "xco2", "xh20", "gden", "wden", "tco2"]
         if dig["case"] != "spe11a":
@@ -708,6 +713,54 @@ def dense_data(dig):
         handle_performance_spatial(dig, dil)
 
 
+def handle_yaxis_mapping(dig, dil):
+    """Extend the indices accounting for the y direction"""
+    simycent = [0.0] * dig["gxyz"][1]
+    with open(f"{dig['path']}/deck/ycenters.txt", "r", encoding="utf8") as file:
+        for j, row in enumerate(csv.reader(file)):
+            simycent[j] = float(row[0])
+    indy = [pd.Series(np.abs(dil["refycent"] - y_c)).argmin() for y_c in simycent]
+    tmp_inds = np.zeros(dig["nocellst"], dtype=int)
+    mults = np.zeros(dig["gxyz"][0], dtype=int)
+    for indz in range(dig["gxyz"][2]):
+        i_i = dig["gxyz"][0] * (dig["gxyz"][2] - indz - 1)
+        i_f = dig["gxyz"][0] * (dig["gxyz"][2] - indz)
+        iii = dig["gxyz"][0] * dig["gxyz"][1] * (dig["gxyz"][2] - 1 - indz)
+        if indz != 0:
+            mults += 1 * (
+                dil["cell_ind"][i_f : i_f + dig["gxyz"][0]]
+                != dil["cell_ind"][i_i : i_i + dig["gxyz"][0]]
+            )
+        values = dil["cell_ind"][i_i : i_i + dig["gxyz"][0]] + mults * (
+            dig["nxyz"][0]
+        ) * (dig["nxyz"][1] - 1)
+        tmp_inds[iii : iii + dig["gxyz"][0]] = values
+        for j, iy in enumerate(indy[1:]):
+            tmp_inds[
+                iii + dig["gxyz"][0] * (j + 1) : iii + dig["gxyz"][0] * (j + 2)
+            ] = (iy * dig["nxyz"][0]) + values
+    dil["cell_ind"] = tmp_inds
+    return dil
+
+
+def handle_inactive_mapping(dig, dil):
+    """Set to inf the inactive grid centers in the reporting grid"""
+    var_array = np.empty(dig["noxz"]) * np.nan
+    var_array[dig["actind"]] = 0.0
+    for i in np.unique(dil["cell_ind"]):
+        inds = i == dil["cell_ind"]
+        if np.isnan(np.sum(var_array[inds])):
+            dil["refxgrid"][i] = np.inf
+            dil["refzgrid"][i] = np.inf
+    ind, dil["cell_ind"] = 0, np.zeros(dig["nocellst"], dtype=int)
+    for i, k in zip(dil["simxcent"], dil["simzcent"]):
+        dil["cell_ind"][ind] = pd.Series(
+            np.abs(dil["refxgrid"] - i) + np.abs(dil["refzgrid"] - k)
+        ).argmin()
+        ind += 1
+    return dil
+
+
 def handle_performance_spatial(dig, dil):
     """Create the performance spatial maps"""
     dil = static_map_to_report_grid_performance_spatial(dig, dil)
@@ -924,27 +977,6 @@ def compute_xh20(dig, dil, h2o_v, co2_g):
     return dil
 
 
-def handle_inactive_mapping(dig, dil):
-    """Set to inf the inactive grid centers in the reporting grid"""
-    var_array = np.empty(dig["nocellst"]) * np.nan
-    var_array[dig["actind"]] = 0.0
-    for i in np.unique(dil["cell_ind"]):
-        inds = i == dil["cell_ind"]
-        if np.isnan(np.sum(var_array[inds])):
-            dil["refxgrid"][i] = np.inf
-            dil["refygrid"][i] = np.inf
-            dil["refzgrid"][i] = np.inf
-    ind, dil["cell_ind"] = 0, np.zeros(dig["nocellst"], dtype=int)
-    for i, j, k in zip(dil["simxcent"], dil["simycent"], dil["simzcent"]):
-        dil["cell_ind"][ind] = pd.Series(
-            (dil["refxgrid"] - i) ** 2
-            + (dil["refygrid"] - j) ** 2
-            + (dil["refzgrid"] - k) ** 2
-        ).argmin()
-        ind += 1
-    return dil
-
-
 def map_to_report_grid(dig, dil, names):
     """Map the simulation grid to the reporting grid"""
     for name in names: