Relative permeability model #1531

wenqing · 2016-11-09T14:30:42Z

This PR presents:

base class for relative permeability models
class for WettingPhaseVanGenuchten model
class for NonWettingPhaseVanGenuchten model
class for WettingPhaseBrookCoreyOilGas model
class for NonWettingPhaseBrookCoreyOilGas model
class for curve type model
unit tests for these classes.

chleh

Only minor things from my side. In general the code looks good.
For general discussion: We have curves globally in the project data. Should they be exclusively parsed there? Should they be removed from there? Should they be passed to all configuration parsing code?

chleh · 2016-11-09T17:12:37Z

...ib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.cpp

+            "The exponent parameter of WettingPhaseVanGenuchten relative\n"
+            " permeability model, m, must be in an interval of [0, 1]");
+    }
+    //! \ogs_file_param{material_property__porous_medium__porous_medium__relative_permeability__WettingPhaseVanGenuchten__m}


Typo: __m at the end?

Corrected to __krel_min

chleh · 2016-11-09T17:12:56Z

...ib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.cpp

+            " permeability model, m, must be in an interval of [0, 1]");
+    }
+    //! \ogs_file_param{material_property__porous_medium__porous_medium__relative_permeability__WettingPhaseVanGenuchten__m}
+    const double krel_max = config.getConfigParameter<double>("krel_min");


Typo: krel_max vs. min?

Corrected to krel_min

chleh · 2016-11-09T17:13:34Z

...ib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.cpp

+            " permeability model, m, must be in an interval of [0, 1]");
+    }
+    //! \ogs_file_param{material_property__porous_medium__porous_medium__relative_permeability__NonWettingPhaseVanGenuchten__m}
+    const double krel_max = config.getConfigParameter<double>("krel_min");


Typo? max vs min? And __m above.

chleh · 2016-11-09T17:14:28Z

...ib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.cpp

+            "relative permeability model, m, must not be smaller than 1");
+    }
+    //! \ogs_file_param{material_property__porous_medium__porous_medium__relative_permeability__WettingPhaseBrookCoreyOilGas__m}
+    const double krel_max = config.getConfigParameter<double>("krel_min");


Looks like a recurrent pattern :-)

chleh · 2016-11-09T17:17:34Z

...ib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.cpp

+/**
+    \param config ConfigTree object which contains the input data
+                  including <type>NonWettingPhaseBrookCoreyOilGas</type>
+                  and it has a tag of <relative_permeability>


Do you mean with that that the XML looks like that:

<relative_permeability> <type>NonWettingPhaseBrookCoreyOilGas</type> ... <relative_permeability>

?

Yes. The unit test has the syntax.

chleh · 2016-11-09T17:28:59Z

...ialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.cpp

+    const double S = MathLib::limitValueInInterval(
+        saturation_w, _Sr + _minor_offset, _Smax - _minor_offset);
+    const double Se = (S - _Sr) / (_Smax - _Sr);
+    const double krel = std::pow(1.0 - Se, 1.0 / 3.0) *


Minor: there is also std::cbrt(), which is probably a bit faster.

chleh · 2016-11-09T17:29:54Z

MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.h

+
+    /// Get relative permeability value.
+    /// \param saturation_w Non-wetting phase saturation
+    double getValue(const double saturation_w) const override;


Likewise: Are derivatives needed?

chleh · 2016-11-09T17:30:47Z

MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/ReletivePermeabilityCurve.h

+{
+namespace PorousMedium
+{
+class ReletivePermeabilityCurve final : public RelativePermeability


typo: Relative

chleh · 2016-11-09T17:32:40Z

MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseVanGenuchten.cpp

+    const double Se = (S - _Sr) / (_Smax - _Sr);
+    const double krel =
+        std::sqrt(Se) *
+        std::pow(1.0 - std::pow(1.0 - std::pow(Se, 1.0 / _mm), _mm), 2);


Minor: std::pow for squaring is rather slow.

chleh · 2016-11-09T17:34:32Z

Tests/MaterialLib/TestRelativePermeabilityModel.cpp

+                             BaseLib::ConfigTree::onwarning);
+    auto const& sub_config = conf.getConfigSubtree("relative_permeability");
+    //! \ogs_file_attr{prj__material_property__porous_medium__porous_medium__id}
+    auto const id = sub_config.getConfigAttributeOptional<int>("id");


If you don't need an attribute, you can use sub_config.ignoreAttribute("id"). I think such a method exists.
Btw: What is the id for? You don't read it in the rest of this PR, so at least you could remove it from the test.

id is reserved for material IDs of elements. We can not input the model parameters as Parameter data because that we may use different models for different material zones in a modelling.

chleh · 2016-11-10T15:20:22Z

Thanks for the update.

wenqing · 2016-11-10T15:31:14Z

@chleh For the curve, when I started to implement the curve type material model, I found that the curve should be placed to where it used or belong to. Parsing all curves and saved them as a member of ProjectData is not the best choice. Material models can have curves, BCs can have curves. These curves should be encapsulated to the classes that use them.

Besides, input the variables and values in separated vector is not good to perform a manual input.

Therefore, I did not add the parsing of the curve for relative permeability to ProjectData class.

endJunction · 2016-11-10T18:27:20Z

...ib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.cpp

+            curve_config.getConfigSubtreeList("point"))
+        {
+            const auto& point =
+                //! \ogs_file_param{material_property__porous_medium__porous_medium__relative_permeability__Curve_curve__points__data}


I guess the ending should be '__point__data', one 's' too much.

endJunction · 2016-11-10T18:28:26Z

Tests/MaterialLib/TestRelativePermeabilityModel.cpp

+        "           <point><data> 0. 0.9   </data></point>"
+        "           <point><data> 0.4 0.5  </data></point>"
+        "           <point><data> 0.9 0.01 </data></point>"
+        "       </curve>"


Why not use the already existing parser for curves? Why do we need another format?

See my message to Christoph:

For the curve, when I started to implement the curve type material model, I found that the curve should be placed to where it used or belong to. Parsing all curves and saved them as a member of ProjectData is not the best choice. Material models can have curves, BCs can have curves. These curves should be encapsulated to the classes that use them. Besides, input the variables and values in separated vector is not good to perform a manual input. Therefore, I did not add the parsing of the curve for relative permeability to ProjectData class.

If you and others think it is not good to have this XML syntax for a curve. I can use the two array type input by extracting createPiecewiseLinearInterpolation(BaseLib::ConfigTree const& config) from Project Data and saving it to a single file in MathLib/InterpolationAlgorithms.

It would be good to have same syntax for the input of similar (if not equal) data. The project files are not straight forward to understand already, but confusing others with two different options for data input (for curves) is not goining to be an improvement.
I don't mind which syntax is used but it should be consistent.

Personally I prefer the coordinates-values input in long lines, because the splitting of the data into coordinate-value pairs does not carry any additional information.

I also would prefer one implementation for reading curves.

@TomFischer I prepared one in PR #1529.
@endJunction @TomFischer. If all of you prefer two array to read a curve. I can change it in #1529.

For this PR. I think, I can drop the curve type and move it a new PR after the curve input is unified.

endJunction · 2016-11-10T18:29:23Z

Tests/MaterialLib/TestRelativePermeabilityModel.cpp

+        if (i == S.size() - 1)
+            continue;
+        // Compare the derivative with numerical one.
+        ASSERT_NEAR((perm_model->getValue(S[i] + purterbation) -


typo in perturbation.

endJunction · 2016-11-10T18:30:48Z

Tests/MaterialLib/TestRelativePermeabilityModel.cpp

+                             BaseLib::ConfigTree::onwarning);
+    auto const& sub_config = conf.getConfigSubtree("relative_permeability");
+    sub_config.ignoreConfigAttribute("id");
+    return MaterialLib::PorousMedium::createRelativePermeabilityModel(


If you are 'using namespace MaterialLib::PorousMedium' than this can be shorter. Choose one variant.

endJunction · 2016-11-10T18:33:04Z

MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.h

+namespace PorousMedium
+{
+/**
+ *   \brief BrookCorey oil-gas model: non-wetting phase


This one seems to be of some van Genuchten type.

endJunction · 2016-11-10T18:34:24Z

...ialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.cpp

+    return (-std::pow(temp_val, 2. * _mm) / (3. * cbrt1_Se * cbrt1_Se) -
+            2. * cbrt1_Se * std::pow(temp_val, 2. * _mm - 1.) *
+                std::pow(Se, (1. - _mm) / _mm)) /
+           (_Smax - _Sr);


Unfortunately (and I made that mistake recently too) identifiers starting with an underscore followed by a capital case letter are reserved for the compilers.

Thanks. Just checked and found names of __ and _[Capital letter] are reserved. I remembered that I've ever changed one member name with _[Capital letter] in MathLib for compilation error but did not check the reason.

I will replace the similar names in the capillary pressure classes in a PR following up PR #1529.

TomFischer · 2016-11-11T11:05:36Z

...ib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.cpp

+
+    //! \ogs_file_param{material_property__porous_medium__porous_medium__relative_permeability__WettingPhaseVanGenuchten__m}
+    const double m = config.getConfigParameter<double>("m");
+    if (m > 1.0)  // m <= 1


If m is negative is not checked.

Added m<0 in the if-condition

TomFischer · 2016-11-11T11:50:04Z

...ib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.cpp

+                  including <type>NonWettingPhaseVanGenuchten</type>
+                  and it has a tag of <relative_permeability>
+*/
+std::unique_ptr<RelativePermeability> createNonWettingPhaseVanGenuchten(


createNonWettingPhaseVanGenuchten and createWettingPhaseVanGenuchten seems to be almost identical. Is it possible to refactor common things?
Update Also the following two create functions seems to be almost identical. All functions read the parameters Sr, Smax, m, and krel_min. The rading and range checking can be unified.

At the beginning, I had only one local function with an option for creating different models. Because of Dima's comment in PR #1517, which has a similar function, I split the single function to four functions. Another consideration is keyword documentation. Although parameter names are the same in the four models, they should be documented for each model individually.

TomFischer · 2016-11-11T12:08:23Z

...rialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/WettingPhaseBrookCoreyOilGas.h

+    const double _saturation_max;  ///< Maximum saturation.
+    const double _mm;              ///< Exponent (>=1.0), n=1/(1-mm).
+    const double _krel_min;        ///< Minimum relative permeability
+};


The above attributes are use in other classes, too (for instance in NonWettingPhaseBrookCoreyOilGas). Maybe it makes sense that these classes have a common base class?

Additional, it would be possible to implement limitValueInInterval(saturation) as a method in this base class. The base class method could use the needed attributes (_saturation_r, _minor_offset, and _saturation_max) and the code would be more readable in my opinion.

@TomFischer limitValueInInterval is also used in capillary pressure classes. Therefore, it was added to MathTool

Moved two common members, _saturation_r and saturation_max to the base class.

endJunction

When the curves input data question is solved in this or in another but consistent way ⏩

TomFischer · 2016-11-17T10:22:16Z

...ib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.cpp

+
+    //! \ogs_file_param{material_property__porous_medium__porous_medium__relative_permeability__NonWettingPhaseBrookCoreyOilGas__m}
+    const double m = config.getConfigParameter<double>("m");
+    if (m < 1.0)  // m >= 1


In this modell m < 0 is allowed?

No. m<0 is included in m<1.0.

In other models you check if m is in [0,1]. Here you do not do this. So is it allowed to give negative values for m in this model?

m is a double paramerer. If a negative value is given for it, the value is smaller than 1.0 too, and the program gives a fatal error as well. For Brook-Corey models, m>1.

For van Genuchten model m in[0,1].

Now I understand. Thanks!

TomFischer · 2016-11-21T05:50:26Z

When rebase is fine, tests are ok: 👍

wenqing · 2016-11-21T16:00:10Z

Thanks for your comments. If there is not further comment, I am going to merge it.

endJunction

Mostly comments and some style suggestions.
⏩ when you think there is already enough improvement.

endJunction · 2016-11-21T19:01:23Z

...ib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.cpp

+
+    //! \ogs_file_param{material__porous_medium__relative_permeability__WettingPhaseVanGenuchten__m}
+    const double m = config.getConfigParameter<double>("m");
+    if (m < 0. || m > 1.0)  // m <= 1


the comment does not match the expression. I'd remove the comments; the reason is already given in the error message.

removed it here and one another in other line.

endJunction · 2016-11-21T19:27:39Z

MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.h

+    NonWettingPhaseVanGenuchten(const double Snr, const double Snmax,
+                                const double m, const double krel_min)
+        : RelativePermeability(1. - Snmax, 1. - Snr),
+          _mm(m),


just curious, why not _m? It would be consistent with the parameter name.

Changed to _m.

endJunction · 2016-11-21T19:28:23Z

MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/NonWettingPhaseVanGenuchten.h

+    double getdValue(const double saturation_w) const override;
+
+private:
+    const double _mm;        ///< Exponent (<=1.0), n=1/(1-mm).


In the corresponding create function you check if this value is in [0,1] interval; the comment here is different. I can not say which one is correct, but if the checks and comments are kept together it is easier to spot mistakes.

Changed the comment to
///< Exponent m, m in [0, 1], n=1/(1-m).

endJunction · 2016-11-21T19:29:50Z

MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h

+    virtual double getdValue(const double saturation) const = 0;
+
+protected:
+    /** A small number for an offset to set the bound of S, the saturation, such


Sometimes doxygen comments are in /** */ style, sometimes in ///. Choose one, stick to it.

I use /// for comments be less than two lines, and /** */ for those more than two lines.

:) the rule is fine, but this particular one has only two lines ;)
But keep it like it is, I don't mind it too much.

endJunction · 2016-11-21T19:34:16Z

Tests/MaterialLib/TestRelativePermeabilityModel.cpp

+{
+    const char xml[] =
+        // Should be <relative_permeability id="0">
+        "<relative_permeability id=\"0\">"


The comment's statement is true. Probably the comment is no longer needed.

… and corrected the documenations for some keywords.

…ments

with respect to saturation.

according to the comments by Dima.

ogsbot · 2020-06-19T12:37:47Z

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

chleh approved these changes Nov 9, 2016

View reviewed changes

wenqing force-pushed the krel branch from c518e02 to 8a0aced Compare November 10, 2016 15:08

endJunction reviewed Nov 10, 2016

View reviewed changes

TomFischer reviewed Nov 11, 2016

View reviewed changes

wenqing mentioned this pull request Nov 11, 2016

[Alg] Added a class for piece wise linear curve #1529

Merged

endJunction approved these changes Nov 12, 2016

View reviewed changes

endJunction force-pushed the master branch from 6be30e3 to 95d89c8 Compare November 12, 2016 16:25

TomFischer reviewed Nov 17, 2016

View reviewed changes

wenqing force-pushed the krel branch from 5b033a2 to 45db30f Compare November 18, 2016 16:30

wenqing force-pushed the krel branch 2 times, most recently from db2cc2f to 1272e53 Compare November 21, 2016 15:36

endJunction approved these changes Nov 21, 2016

View reviewed changes

endJunction added the please rebase label Nov 21, 2016

wenqing added 6 commits November 22, 2016 10:06

[Unsat] Added the base class for relative permeability

cdcb89c

[Unsat] Added Brook-Corey oil-gas relative permeability model

db90212

[Unsat] Added van Genuchten relative pemeability model

35232b5

[Unsat] Added curve type relative permeability model

f02a133

[Unsat] Added the creator for relative permeability models

06347c8

[Unsat] Add unit test for relative permeability models

5207e3f

wenqing added 9 commits November 22, 2016 10:06

[Unsat] Removed copy constructor and moved a member to base class.

f25d8cc

[Unsat] Changed the parameters of the constructors of several classes…

92b3fa1

… and corrected the documenations for some keywords.

[Unsat] Some corrections and changes according to the Christoph's com…

30c7705

…ments

[Unsat] Added the calculation of the derivative of relative permeability

2d9cff7

with respect to saturation.

[Unsat] Changed some member names and made other minor changes

0261038

according to the comments by Dima.

[Unsat] Check m<0 for van Genuchten relative permeability model

a53bcb1

[Unsat] Moved the two common members to the base class

94afa3e

[Unsat] Changed the XML syntax for curve

bbadb6a

[Unsat] Changed the documentation for relative permeability keywords

6f13cd3

wenqing force-pushed the krel branch from 1336632 to 9734e16 Compare November 22, 2016 09:07

[Unsat] Changed one class member name and some comments

5bcbc36

wenqing force-pushed the krel branch from 9734e16 to 5bcbc36 Compare November 22, 2016 11:45

wenqing merged commit e692786 into ufz:master Nov 22, 2016

endJunction removed the please rebase label Nov 22, 2016

Relative permeability model #1531

Relative permeability model #1531

Conversation

wenqing commented Nov 9, 2016

chleh left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

chleh commented Nov 10, 2016

wenqing commented Nov 10, 2016

Choose a reason for hiding this comment

endJunction Nov 10, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

TomFischer Nov 11, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

endJunction left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

TomFischer commented Nov 21, 2016

wenqing commented Nov 21, 2016

endJunction left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ogsbot commented Jun 19, 2020

endJunction Nov 10, 2016 •

edited

Loading

TomFischer Nov 11, 2016 •

edited

Loading