Fixups for check_swatinit after user testing (#343)

* Fixed palette for scatter plots
* Regenerated and almost reproducible plots in doc
* Adjust placement of annotations in volplot
* Reorder arguments in parser, for --help output
* Allow EQLNUM=1 to not exist without crash
* Ensure large dataframes are printed in full to stdout.
* Add interactive plots to pytest, not run by default
* Add SWU to output csv and notify that SWU<1 is not supported
* Improve docs
* Avoid PC computation when PPCW is dummy (SWAT=SWL)
* Reorder CSV output to prioritize non-nan data

docs/scripts/check_swatinit.rst

@@ -71,13 +71,15 @@ are:
   significant volumes, revise the modelling.
 
 ``SWATINIT_1``
-  When SWATINIT is 1 above the contact, Eclipse will ignore SWATINIT in the cell
-  and not touch the capillary pressure function. This will typically result in extra
-  hydrocarbons added to the model for a normal capillary pressure function. This
-  could be ok as long as the porosities and/or permeabilities of these cells are
-  small. If it is not, you should look into if there are upscaling issues for
-  this cell. In situations with nonzero item #9 in EQUIL, this can also occur
-  below the contact.
+  When SWATINIT is 1 above the contact, Eclipse will ignore SWATINIT in the
+  cell and not touch the capillary pressure function. This will typically
+  result in extra hydrocarbons added to the model for a normal capillary
+  pressure function. This could be ok as long as the porosities and/or
+  permeabilities of these cells are small. If it is not, you should look into
+  if there are upscaling issues for this cell. In situations with nonzero item
+  #9 in EQUIL, this can also occur below the contact. If SWU is included in the
+  model and less than 1, cells where SWATINIT is equal or larger than SWU will
+  also be flagged as SWATINIT_1 as the behaviour is the same.
 
 ``HC_BELOW_FWL``
   If SWATINIT is less than 1 below the contact provided in EQUIL, Eclipse will
@@ -96,12 +98,24 @@ are:
   SWATINIT was 1 in the water zone, and SWAT is set to 1.
 
 
-Example text output
--------------------
+Text output from check_swatinit
+-------------------------------
+
+When run interactively in a terminal on an Eclipse case, check_swatinit will
+print a text output summarizing the results from flagging each cell into
+distinct QC categories.  The table below sums the water volume changes from
+SWATINIT to SWAT for each QC category, in terms of *reservoir* volumes, i.e. no
+conversion to surface conditions. The percentages are with respect the
+SWATINIT_WVOL (i.e. the water volume that SWATINIT requested).
+
+To the right in the output, the corresponding hydrocarbon porevolumes under
+*reservoir* conditions are given. The number 66.571 is SWATINIT_WVOL
+substracted from PORV, and then the percentages below are the volumetric change
+in water volume with respect to this hydrocarbon pore volume.
 
 .. code-block:: console
 
-  $ check_swatinit DROGON.DATA
+  $ check_swatinit DROGON-EXAMPLE.DATA
   VOLUME                     3203.1103 Mrm3
   PORV                        571.1770 Mrm3
   SWATINIT_WVOL               504.6057 Mrm3           HC:   66.571 Mrm3
@@ -123,6 +137,23 @@ Example text output
   1      1       22.626438    0.136527
 
 
+The maximal values tables then gives the maximum capillary pressure as present
+in the input SWOF/SWFN table (the number PCOW_MAX) for each SATNUM.  The number
+PPCW in the next table is the corresponding number after Eclipse has scaled the
+capillary pressure, and you can find here that PPCW = PCOW_MAX * PC_SCALING. In
+the ideal case (no error from upscaling etc), PC_SCALING would be 1 in all
+cells. Since PC_SCALING is less than 1 it means that capillary pressure had to
+be downscaled by at least a factor ~10 in every cell, meaning that the
+capillary pressure input curve should be revisited.
+
+Use also the scatter plot through the ``--plot`` option to get an impression
+of which capillary pressures ranges are present in your reservoir model.
+
+More often, PC_SCALING and PPCW will be large, with scaling larger than 1.
+Still, this table only prints the maximal values, and if severe PC_SCALING only
+occurs in a few cells, it could still be acceptable.
+
+
 Outputted CSV file
 ------------------
 
@@ -159,15 +190,15 @@ Example plots
    :width: 70%
 
    A waterfall chart illustrating what contributes to the change from SWATINIT
-   to SWAT. This plot is obtained by adding the ``--volplot`` command line option.
-   The numbers inside the plot is the percentage change in terms of reservoir
-   volumes, blue numbers are with respect to SWATINIT_WVOL and green numbers
-   are with respect to initial hydrocarbon volumes.
+   to SWAT. This plot is obtained by adding the ``--volplot`` command line
+   option.  The numbers inside the plot are the percentage changes in terms of
+   reservoir volumes, blue numbers are with respect to SWATINIT_WVOL and green
+   numbers are with respect to initial hydrocarbon volumes.
 
 
 .. figure:: images/check_swatinit_scatter.png
    :align: center
-   :width: 90%
+   :width: 98%
 
    A panel of reservoir properties versus depth, coloured by the QC_FLAG, for a
    specific EQLNUM. Use the command line option ``--plot`` together with

docs/scripts/images/make_check_swatinit_images.sh

@@ -0,0 +1,10 @@
+#!/bin/bash
+
+cd "$(dirname "$0")"
+
+check_swatinit ../../../tests/data/reek/eclipse/model/2_R001_REEK-0.DATA \
+    --volplotfile check_swatinit_volplot.png
+
+# Sorry, this csv-file is not under version control:
+check_swatinit ~/drogon.csv --plotfile check_swatinit_scatter.png  \
+    --eqlnum 4

setup.cfg

@@ -26,6 +26,7 @@ addopts =
 
 markers =
     integration: marks a test as an integration test
+    plot: marks a test as interactive, plots will flash to the screen
 
 [tool:pylint]
 # Module docstrings are not required, there are other means of documenting at

src/subscript/check_swatinit/check_swatinit.py

@@ -69,7 +69,7 @@ def main():
 
         if args.output != "":
             logger.info("Exporting CSV to %s", args.output)
-            qc_frame.to_csv(args.output, index=False)
+            reorder_dframe_for_nonnans(qc_frame).to_csv(args.output, index=False)
 
     if "SWATINIT" not in qc_frame:
         print("Model did not use SWATINIT")
@@ -84,7 +84,15 @@ def main():
     print()
     print(human_report_pc_scaling(qc_frame))
 
-    if args.eqlnum not in qc_frame["EQLNUM"].values:
+    if args.volplot or args.volplotfile:
+        plotter.wvol_waterfall(qc_vols)
+    if args.volplot:
+        pyplot.show()
+    if args.volplotfile:
+        print(f"Dumping volume plot to {args.volplotfile}")
+        pyplot.savefig(args.volplotfile)
+
+    if (args.plotfile or args.plot) and args.eqlnum not in qc_frame["EQLNUM"].values:
         sys.exit(f"Error: EQLNUM {args.eqlnum} does not exist in grid. No plotting.")
     if args.plot or args.plotfile:
         plotter.plot_qc_panels(qc_frame[qc_frame["EQLNUM"] == args.eqlnum])
@@ -94,14 +102,6 @@ def main():
         print(f"Dumping plot to {args.plotfile}")
         pyplot.savefig(args.plotfile)
 
-    if args.volplot or args.volplotfile:
-        plotter.wvol_waterfall(qc_vols)
-    if args.volplot:
-        pyplot.show()
-    if args.volplotfile:
-        print(f"Dumping volume plot to {args.volplotfile}")
-        pyplot.savefig(args.volplotfile)
-
 
 def check_applicability(eclfiles):
     """Check that the input is relevant for usage with check_swatinit. This
@@ -140,6 +140,17 @@ def check_applicability(eclfiles):
         )
 
 
+def reorder_dframe_for_nonnans(dframe):
+    """Reorder a dataframe so that rows with less NaN comes first, this
+    will aid data analysis application to deduce correct datatypes for
+    columns"""
+    null_count = "__NULL_COUNT__"
+    dframe[null_count] = dframe.isnull().sum(axis=1)
+    return (
+        dframe.sort_values(null_count).drop(null_count, axis=1).reset_index(drop=True)
+    )
+
+
 def human_report_qc_vols(qc_vols):
     """Produce a string with a human report for volumes
 
@@ -195,12 +206,14 @@ def human_report_pc_scaling(qc_frame):
     string = ""
     string += "Maximal values:\n"
     string += "---------------\n"
-    string += str(qc_frame.groupby("SATNUM").max()[["PCOW_MAX"]])
+    string += qc_frame.groupby("SATNUM").max()[["PCOW_MAX"]].to_string()
     string += "\n"
-    string += str(qc_frame.groupby(["EQLNUM", "SATNUM"]).max()[["PPCW", "PC_SCALING"]])
+    string += (
+        qc_frame.groupby(["EQLNUM", "SATNUM"]).max()[["PPCW", "PC_SCALING"]].to_string()
+    )
     string += "\n\n"
-    string += "EQUIL initialization option #9):\n"
-    string += str(qc_frame.groupby("EQLNUM").max()[["OIP_INIT"]].astype(int))
+    string += "EQUIL initialization option #9:\n"
+    string += qc_frame.groupby("EQLNUM").max()[["OIP_INIT"]].astype(int).to_string()
     return string
 
 
@@ -227,6 +240,7 @@ def make_qc_gridframe(eclfiles):
             "PCW",
             "PPCW",
             "SWL",
+            "SWU",  # It is a feature request to process this
             "SWLPC",  # Extract in case it is there, but it is not supported yet.
         ],
         rstdates="first",
@@ -243,6 +257,9 @@ def make_qc_gridframe(eclfiles):
         logger.warning("Consider adding FILLEPS to the PROPS section")
         grid_df["SWL"] = 0.0
 
+    if "SWU" in grid_df and (grid_df["SWU"] - 1).abs().sum() > 0:
+        logger.warning("SWU is less than 1, not yet supported by check_swatinit")
+
     if "SWLPC" in grid_df and (grid_df["SWL"] - grid_df["SWLPC"]).abs().sum() > 0:
         raise ValueError("SWLPC is in the data, but not supported by check_swatinit")
 
@@ -510,6 +527,14 @@ def compute_pc(qc_frame, satfunc_df):
     else:
         contact = "GWC"
 
+    # When SWATINIT=SWL=SWAT, PPCW as reported by Eclipse is the
+    # same as PCOW_MAX, and we cannot use it to compute PC, remove it:
+    if "SWL" in qc_frame:
+        p_cap[
+            np.isclose(qc_frame["SWAT"], qc_frame["SWL"])
+            & np.isclose(qc_frame["PC_SCALING"], 1)
+        ] = np.nan
+
     if "QC_FLAG" in qc_frame:
         p_cap[
             (qc_frame["QC_FLAG"] == __SWATINIT_1__)
@@ -562,11 +587,7 @@ def merge_equil(grid_df, equil_df):
     assert (
         not pd.isnull(equil_df).any().any()
     ), f"BUG: NaNs in equil dataframe:\n{equil_df}"
-    grid_df = grid_df.merge(
-        equil_df,
-        on="EQLNUM",
-        how="left",
-    )
+    grid_df = grid_df.merge(equil_df, on="EQLNUM", how="left")
     return grid_df
 
 
@@ -626,6 +647,19 @@ def get_parser():
         default="",
         help="Output filename for CSV that can be used for QC in other tools",
     )
+    parser.add_argument(
+        "--volplot",
+        action="store_true",
+        help=(
+            "Display a waterfall chart with the water "
+            "saturation volumes from SWATINIT to SWAT"
+        ),
+    )
+    parser.add_argument(
+        "--volplotfile",
+        type=str,
+        help="PNG filename for where to dump a waterfall chart.",
+    )
     parser.add_argument(
         "--plot",
         action="store_true",
@@ -634,6 +668,9 @@ def get_parser():
             "region, use together with the --eqlnum option"
         ),
     )
+    parser.add_argument(
+        "--plotfile", type=str, help="PNG filename for where to dump a QC plot."
+    )
     parser.add_argument(
         "--eqlnum",
         type=int,
@@ -643,22 +680,6 @@ def get_parser():
             "Does not affect CSV output"
         ),
     )
-    parser.add_argument(
-        "--plotfile", type=str, help="PNG filename for where to dump a QC plot."
-    )
-    parser.add_argument(
-        "--volplot",
-        action="store_true",
-        help=(
-            "Display a waterfall chart with the water "
-            "saturation volumes from SWATINIT to SWAT"
-        ),
-    )
-    parser.add_argument(
-        "--volplotfile",
-        type=str,
-        help="PNG filename for where to dump a waterfall chart.",
-    )
     return parser