Latent issue in EndpointInitializer::findVerticalPoints()

[bska]

Lines 197--212 of SatfuncPropertyInitializers.hpp (method findVerticalPoints() of class EndpointInitializer) contains the following code

for (size_t rowIdx = 0; rowIdx < krwCol.size(); ++rowIdx) {
    if (krwCol[rowIdx] == 0.0) {
        m_krorw[tableIdx] = krowCol[rowIdx - 1];
        break;
    }
}

// find the oil relperm which corresponds to the critical gas saturation
const auto &krgCol = sgofTables[tableIdx].getKrgColumn();
const auto &krogCol = sgofTables[tableIdx].getKrogColumn();
for (size_t rowIdx = 0; rowIdx < krgCol.size(); ++rowIdx) {
    if (krgCol[rowIdx] == 0.0) {
        m_krorg[tableIdx] = krogCol[rowIdx - 1];
        break;
    }
}

The expressions rowIdx - 1 produce out-of-bounds access ((size_t) -1) when rowIdx == 0 (i.e., on the first iteration). I don't know what the correct fix is in this case.

[andlaus]

referring to the Eclipse technical description, the case you're thinking about is most likely incorrect data for the satfunc keywords. that's because the TD says that the critical saturation of a phase is the highest saturation where the corresponding relperm is zero.

this definition is obviously non-sense if there is no such saturation. I agree that this case should be handled, though. exceptions anyone?

[bska]

referring to the Eclipse technical description, the case you're thinking about is most likely incorrect data for the satfunc keywords.

Possibly, but the check is inverted in that case. If the objective is to locate the critical saturation you need to look for the first non-zero relative permeability and only extract the previous row index when you've verified that it actually exists. The way this check is written you are guaranteed to reference table[(size_t) -1] if there is any zero relative permeabilities at the beginning of the table.

[bska]

If the objective is to locate the critical saturation you need to look for the first non-zero relative permeability and only extract the previous row index when you've verified that it actually exists.

Just to follow up on this, here is MRST's table_points function for SWOF. This function computes connate water saturation (Swco), critical water saturation (Swcr), critcal oil saturation in oil/water system (Sowcr), maximum water saturation (Swmax), and relative permeability for oil at connate water (krocw) for a single SWOF table T represented as an m-by-4 double array. MRST also provides similar functions for keywords SGOF, SWFN, SGFN, SOF2, and SOF3.

function [Swco, Swcr, Sowcr, Swmax, krocw] = table_points(T)
   % Basic validation of input data.
   assert (~any(any(T(:, 1:3) < 0)), ...
           'Saturation and rel-perm values must be non-negative.');

   assert (~ (T( 1 ,2) > 0), ...
           'Water must be immobile at connate water saturation.');
   assert (~ (T(end,3) > 0), ...
           'Oil must be immobile at maximum water saturation.');

   assert (all (diff(T(:,1)) > 0), ...
           'Water saturation must be monotonically increasing.');
   assert (~any(diff(T(:,2)) < 0), ...
           'Water rel-perm must be level or increasing down column.');
   assert (~any(diff(T(:,3)) > 0), ...
           'Oil rel-perm must be level or decreasing down column.');

   %-----------------------------------------------------------------------
   % Data OK.  Extract critical values from table.
   %

   % 1) Connate water saturation (>= 0) is first Sw encountered in table.
   %
   Swco = T(1,1);

   % 2) Critical water saturation (>= 0) is highest Sw for which krw(Sw)=0.
   %
   i = find(T(:,2) > 0, 1, 'first');
   assert (~isempty(i), ...
           'Table must define a critical water saturation, Swcr.');

   Swcr = T(i - 1, 1);

   % 3) Critical oil saturation (>= 0) is highest So for which kro(So)=0.
   %
   %      Sowcr = max { S_o | k_{row}(S_o) = 0 },  S_o \in [0,1]
   %
   %    This is 1-lowest(Sw) at which the same criterion is satisfied.
   %
   i     = find(~(T(:,3) > 0), 1, 'first');
   Sowcr = 1 - T(i, 1);

   % 4) Maximal water saturation in table.
   Swmax = T(end, 1);

   % 5) Oil relperm at Sw = Swco (=> So = 1 - Swco)
   %
   krocw = T(1, 3);
end

[andlaus]

you are guaranteed to reference table[(size_t) -1] if there is any zero relative permeabilities at the beginning of the table.

True. the point is that a non-zero relperm for a saturation of zero will blow the simulation up in more subtle ways anyway -- if for nothing else, then because it will allow fluxes of a phase which does not exist at a given time...

I'd propose to fix it this way:

if (rowIdx == 0 || kr[rowIdx - 1] != 0)
  throw std::runtime_error("Invalid data specified for keyword.");

[bska]

the point is that a non-zero relperm for a saturation of zero will blow the simulation up in more subtle ways anyway

We're apparently talking past each other. Consider the SWOF table from opm-core's satfuncEPS_A.DATA

SWOF 
0.1   0.0   1.0   0.9
0.2   0.0   0.8   0.8
0.3   0.1   0.6   0.7
0.4   0.2   0.4   0.6
0.7   0.5   0.1   0.3
0.8   0.6   0.0   0.2
0.9   0.7   0.0   0.1
/

This is a perfectly valid saturation function which says that kr_w = 0 when S_w <= 0.2. Consequently the critical water saturation is 0.2. However, the code ends up referencing krw[(size_t) -1] because krw[0] = 0. That's a bug, because krw[0] = 0 is supposed to be there in well-formed input. In fact opm-core's test_satfunc fails due to index out of bounds when using libstdc++'s debug mode.

In short, your search algorithm is wrong. Starting at index one rather than zero may be sufficient, but I'd rather you employ a different strategy altogether if the objective is to compute the critical water saturation.

[andlaus]

That's a bug

Oh, you're right. the condition in the ifs should be the other way around:

if (krwCol[rowIdx] > 0.0)
  // ...

it is not a latent bug though. Will do a PR later today...

[alfbr]

throw std::runtime_error("Invalid data specified for keyword.");

Please no, with lots of cream and sugar on top. While a simulator may want to throw on that input (if it simulates multiphase flow), a 3D viewer most certainly do not. I suggest this is handled by sending an error to the logger.