OPM / opm-parser

http://www.opm-project.org
11 stars 44 forks source link

Bug in 3D-property initialization #1017

Open andlaus opened 7 years ago

andlaus commented 7 years ago

edit: although there seems to be a bug related to the property initialization, the initial analysis is not correct. see below: https://github.com/OPM/opm-parser/issues/1017#issuecomment-269038613

The Norne deck from opm-data features the "SWL" and "ISWL" double grid properties and the latter is supposed to be identical to the former because it is copied via the COPY statement on line 342 of the deck's main file. Unfortunately opm-parser currently thinks otherwise: Running the following program on the Norne deck

#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>

#include <iostream>
#include <tuple>
#include <vector>

int main(int argc, char** argv)
{
    if (argc != 2) {
        std::cout << "Usage: " << argv[0] << " FILNAME\n";
        return 1;
    }

    Opm::Parser parser;
    typedef std::pair<std::string, Opm::InputError::Action> ParseModePair;
    typedef std::vector<ParseModePair> ParseModePairs;
    ParseModePairs tmp;
    tmp.push_back(ParseModePair(Opm::ParseContext::PARSE_RANDOM_SLASH , Opm::InputError::IGNORE));
    Opm::ParseContext parseContext(tmp);

    const Opm::Deck deck = parser.parseFile(argv[1], parseContext);
    const Opm::EclipseState eclState(deck, parseContext);

    if (!eclState.get3DProperties().hasDeckDoubleGridProperty("SWL"))
        std::cout << "PARSER BUG: the Norne deck features the SWL property!\n";

    if (!eclState.get3DProperties().hasDeckDoubleGridProperty("ISWL"))
        std::cout << "PARSER BUG: the Norne deck features the ISWL property!\n";

    const std::vector<double> &swl =
        eclState.get3DProperties().getDoubleGridProperty("SWL").getData();

    const std::vector<double> &iswl =
        eclState.get3DProperties().getDoubleGridProperty("ISWL").getData();

    for (unsigned i = 0; i < swl.size(); ++i) {
        if (swl[i] != iswl[i]) {
            std::cout << "PARSER BUG: In the Norne deck the SWL and ISWL properties are identical!\n";
            std::cout << "(For cell " << i << " opm-parser says SWL is " << swl[i] << " and ISWL is " << iswl[i] << ".)\n";
            break;
        }
    }

    return 0;
}

yields

PARSER BUG: In the Norne deck the SWL and ISWL properties are identical!
(For cell 0 opm-parser says SWL is 0.04 and ISWL is 0.)
joakim-hove commented 7 years ago

Well - there might certainly be bugs, but I do not agree with the current analysis. The relevant part of the COPY statement is:

COPY
 ...
 'SWL'   'ISWL'    1 46 1 112  5  22 /
 ...
/

i.e. only the layers k ∈ [4,21] are copied. So - the difference between SWL and ISWL for the first elements is legitimate. Since the ISWL is not explicitly assigned to in the deck it will assume default values for the layers k ∈ [0,3].

joakim-hove commented 7 years ago

As I wrote in my previous comment I do not think this is a bug; but of course please elaborate and give more details if you think I misunderstand the case.

andlaus commented 7 years ago

I suppose you're right. That said, I'm pretty sure this was not intended by the deck's authors: it makes the endpoints of imbibition and drainage for hysteresis mismatch which is disallowed by the documentation as far as I can remember.

anyway, closing.

joakim-hove commented 7 years ago

I'm pretty sure this was not intended by the deck's authors:

That I don't know of course; but there are no active cells in the top four layers so my hunch is that the original deck authors used some ill-advised assignment saving algorithm.

andlaus commented 7 years ago

but there are no active cells in the top four layers

really? the reason I've stumbled over this was that the material laws initialize differently for the second active cell. maybe I have another bug, though...

joakim-hove commented 7 years ago

really?

Ehhh no - sorry; my bad. There are certainly active cells in the four top layers.

When looking at the INIT file from ECLIPSE the first values of the ISWL keyword are -100000002004087734272.000 - that could be the sign of some "magic" unitialized value which eclipse makes sensible use of? So - there might still be a bug (or maybe more precisely: not-correctly- implemented-undefined-behavior-in-the-case-of-invalid-user-input)? Feel free to reopen if you feel appropriate - although I think the initial analysis was to simple.

image

andlaus commented 7 years ago

Feel free to reopen if you feel appropriate - although I think the initial analysis was to simple.

ok, I'll reopen it as a reminder. Merry Christmas, too :)

OPMUSER commented 7 years ago

Copy should over write layers 5 to 22, why do you say 4 to 21? ThunderBird Signature File

Regards, David Baxendale Email: David.Baxendale@gmail.com


On 2016-12-24 03:17, Joakim Hove wrote:

Well - there might certainly be bugs, but I do not agree with the current analysis. The relevant part of the |COPY| statement is:

|COPY ... 'SWL' 'ISWL' 1 46 1 112 5 22 / ... / |

i.e. only the layers k ∈ [4,21] are copied. So - the difference between |SWL| and |ISWL| for the first elements is legitimate. Since the |ISWL| is not explicitly assigned to in the deck it will assume default values for the layers k ∈ [0,3].

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/OPM/opm-parser/issues/1017#issuecomment-269033914, or mute the thread https://github.com/notifications/unsubscribe-auth/ARIID6mNWKyJrpNfOWTECgt4EgCXYkhKks5rLB4vgaJpZM4LU5V0.

joakim-hove commented 7 years ago

OPM. uses zero offset variables internally - I just did that subtraction upfront.

OPMUSER commented 7 years ago

Joakim,

I'm confused ,the mods should apply to layers 5 to 22 in the user domain, is that equivalent to 4 to 21 in OPM domain? If this was tranmissibility then it would either be 4 to 21 or 6 to 23 depending on the numerical scheme.

jokva commented 7 years ago

Internally, OPM represents layer 1 with numerical zero.

andlaus commented 7 years ago

The confusion comes from FORTRAN vs C style indices: FORTRAN (and the ECL file format) uses the indices 1, .., n to access arrays of size n, in C it is 0, ..., n-1. since OPM is written in C(++) and Eclipse in FORTRAN there's this small "impedance mismatch".

OPMUSER commented 7 years ago

Thanks for the clarification. ThunderBird Signature File

Regards, David Baxendale Email: David.Baxendale@gmail.com


On 2016-12-24 16:44, jokva wrote:

Internally, OPM represents layer 1 with numerical zero.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/OPM/opm-parser/issues/1017#issuecomment-269075510, or mute the thread https://github.com/notifications/unsubscribe-auth/ARIIDwTifMde4bzloehiMFah8m8OvmnMks5rLNtugaJpZM4LU5V0.

OPMUSER commented 7 years ago

Andreas thanks for the clarification, being an old FORTRAN guy, you can see why I was confused :-) ThunderBird Signature File

Regards, David Baxendale Email: David.Baxendale@gmail.com


On 2016-12-24 17:24, Andreas Lauser wrote:

The confusion comes from FORTRAN vs C style indices: FORTRAN (and the ECL file format) uses the indices 1, .., n to access arrays of size n, in C it is 0, ..., n-1. since OPM is written in C(++) and Eclipse in FORTRAN there's this small "impedance mismatch".

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/OPM/opm-parser/issues/1017#issuecomment-269076739, or mute the thread https://github.com/notifications/unsubscribe-auth/ARIID8T6gimgW-7io97RxNBuRpNvujmSks5rLOTVgaJpZM4LU5V0.