Capillary pressure models #1517

wenqing · 2016-11-03T10:44:08Z

This PR presents:

Base class for capillary pressure model
class for van Genuchten capillaryPressure model
class for BrookCoreyCapillaryPressure model
Creator for capillary pressure models
unit tests for these classes

endJunction · 2016-11-03T11:58:50Z

MathLib/MathTools.h

@@ -120,6 +120,16 @@ inline double to_radians(double degrees) {
    return degrees*boost::math::constants::pi<double>()/180.;
 }

+template<typename Type> Type limitValueInInterval(const Type variable,
+                                               const Type lower_bound,
+                                               const Type upper_bound)


~~Better to pass the types by universal reference (Type&&) to avoid copies if the types are expensive to copy.~~
My bad. Not working this way.

The type could be integer or real type.

endJunction · 2016-11-03T12:01:07Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.h

+    double getCapillaryPressure(const double saturation) const override;
+
+    /// Get capillary pressure.
+    double getSturation(const double capillary_pressure) const override;


Comment mentions pressure, but this is saturation.
Typo in saturation.

endJunction · 2016-11-03T12:02:55Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.h

+ *    Note:
+ *     \f[m=1/(1-n)\f]
+ */
+class BrookCoreyCapillaryPressure final : public CapillaryPressureSaturation


Maybe better to call it BrookCoreyCapillaryPressureSaturation to reflect that it is having both p and s and also is derived from CapillaryPressureSaturation.

Changed the name accordingly.

endJunction · 2016-11-03T12:04:53Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.h

+    {
+    }
+
+    BrookCoreyCapillaryPressure(const BrookCoreyCapillaryPressure& orig)


This can be shortened to

BrookCoreyCapillaryPressure(const BrookCoreyCapillaryPressure&) = default;

Also if copy constructor is present, add the move ctor, and maybe both assignments.

Removed the copy constructor because it is not used.

endJunction · 2016-11-03T12:07:39Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressure.h

+    const double _mm;      ///< Exponent (<=1.0), n=1/(1-mm).
+    const double _Pc_max;  ///< Maximum capillaray pressure
+
+    const double _perturbation = std::numeric_limits<double>::epsilon();


Please add a short comment how the _perturbation is used.
From what I have seen in the code, maybe 'offset' would be a better name.

Changed the name to _minor_offset

endJunction

Minor changes required. Looking good otherwise. I didn't check the formulae.

endJunction · 2016-11-03T12:08:19Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h

+    virtual double getCapillaryPressure(const double saturation) const = 0;
+
+    /// Get capillary pressure.
+    virtual double getSturation(const double capillary_ressure) const = 0;


same with comment and typo here.

endJunction · 2016-11-03T12:13:14Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp

+*/
+static std::unique_ptr<CapillaryPressureSaturation>
+createBrookCoreyOrVanGenuchten(BaseLib::ConfigTree const& config,
+                               const bool is_brook_corey)


I'd split it into two different create functions and avoid the boolean parameter.
Adding a third model here would be difficult.

wenqing · 2016-11-03T15:36:47Z

@endJunction Formulas are checked by the unit tests, which use the results from OGS5 as expected data for comparison.

endJunction · 2016-11-03T15:51:02Z

...b/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h

@@ -89,9 +82,9 @@ class BrookCoreyCapillaryPressure final : public CapillaryPressureSaturation
    const double _mm;      ///< Exponent (<=1.0), n=1/(1-mm).
    const double _Pc_max;  ///< Maximum capillaray pressure

-    const double _perturbation = std::numeric_limits<double>::epsilon();
+    const double _minor_offset = std::numeric_limits<double>::epsilon();


Still needs a description how it is used.

Added a description.

endJunction · 2016-11-03T15:52:45Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp

+*/
+static std::unique_ptr<CapillaryPressureSaturation>
+                         createVanGenuchten(BaseLib::ConfigTree const& config,
+                                            const bool is_brook_corey)


boolean is no longer needed, isn't it?

ogsbot · 2016-11-03T17:37:45Z

Jenkins: OGS-6/Linux-PRs-dynamic failed: https://svn.ufz.de:8443/job/OGS-6/job/Linux-PRs-dynamic/1579/

…a's comment

…argument

chleh

Minor comments from my side.

chleh · 2016-11-08T13:12:04Z

...PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.cpp

+    const double saturation) const
+{
+    const double S = MathLib::limitValueInInterval(
+        saturation, _Sr + _minor_offset, _Smax - _minor_offset);


Are there any mathematical concerns here? p_C(S) is cut off at certain values, but the derivative does not take that into account.
Maybe that's common practice, though...

It is used to avoid the overflow in std::pow(Se, -n) when S=0. When S=Smax, setting S=Smax - _minor_offset avoids Pc=0. Therefore, it is a numerical 'tricky' to keep computation run with a tiny losing in accuracy.

Thank you for the explanation.

chleh · 2016-11-08T13:13:39Z

...b/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h

+     *                     [4] \f$ P_c^{\mbox{max}}\f$
+     */
+    BrookCoreyCapillaryPressureSaturation(
+        std::array<double, 5> const& parameters)


I'd prefer to pass the parameters individually rather than as array, because they all mean different things.

chleh · 2016-11-08T13:15:38Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp

+{
+/**
+    \param config ConfigTree object which contains the input data
+                  including <type>BrookCorey</type>


I'd put a check that the type is really BrookCorey.

chleh · 2016-11-08T13:16:02Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp

+
+/**
+    \param config ConfigTree object which contains the input data
+                  including <type>vanGenuchten</type>


Likewise: I'd put a check to the function body.

chleh · 2016-11-08T13:17:03Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp

+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__type}
+    auto const type = config.getConfigParameter<std::string>("type");


Note: If you'd add checks to the specific create functions, you'd switch to peek...() here.

@chleh Thanks for above comments. I will improve the classes in the Fluid directories according to these comments once the current work is finished.

chleh · 2016-11-08T13:18:00Z

...rousMedium/UnsaturatedProperty/CapillaryPressure/vanGenuchtenCapillaryPressureSaturation.cpp

+ *  Created on October 28, 2016, 6:05 PM
+ */
+
+#include "vanGenuchtenCapillaryPressureSaturation.h"


Almost all files start uppercase. Maybe we should stick to this convention as much as possible.

chleh · 2016-11-08T13:19:07Z

...PorousMedium/UnsaturatedProperty/CapillaryPressure/vanGenuchtenCapillaryPressureSaturation.h

+     *                     [4] \f$ P_c^{\mbox{max}}\f$
+     */
+    vanGenuchtenCapillaryPressureSaturation(
+        std::array<double, 5> const& parameters)


As before: I'd prefer no array to be used.

chleh · 2016-11-08T13:21:00Z

Tests/MaterialLib/TestCapillaryPressureSaturationModel.cpp

+    for (std::size_t i = 0; i < S.size(); i++)
+    {
+        ASSERT_NEAR(pc[i], pc_model->getCapillaryPressure(S[i]), 1.e-5);
+        ASSERT_NEAR(S[i], pc_model->getSaturation(pc[i]), 1.e-5);


That doesn't seem too strict, since S is on the order of 0.1...1
Or is it impossible to decrease this tolerance.

Decreased the tolerance for S to 1e-14.

chleh · 2016-11-08T13:21:34Z

Tests/MaterialLib/TestCapillaryPressureSaturationModel.cpp

+    for (std::size_t i = 0; i < S.size(); i++)
+    {
+        ASSERT_NEAR(pc[i], pc_model->getCapillaryPressure(S[i]), 1.e-5);
+        ASSERT_NEAR(S[i], pc_model->getSaturation(pc[i]), 1.e-5);


Maybe it's possible to make the tolerance a bit stricter?

Decreased the tolerance for S to 1e-14.

chleh · 2016-11-08T13:23:02Z

...PorousMedium/UnsaturatedProperty/CapillaryPressure/vanGenuchtenCapillaryPressureSaturation.h

+ * as
+ *    \f[p_b=\rho g/\alpha\f]
+ */
+class vanGenuchtenCapillaryPressureSaturation final


Class names should start uppercase.

added model type check in the function body according to the comments by Christoph

endJunction · 2016-11-08T23:09:04Z

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp

+    //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__type}
+    config.checkConfigParameter("type", "BrookCorey");
+
+    //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__type__pd}


If I'm not mistaken the type in the string must be BrookCorey because it is specific for that model.
@chleh Please correct me if I'm talking garbage.

You're right. The line must be
//! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__BrookCorey__pd}
same for the other doxygen comments.
IMHO since documenting input file parameters is somewhat tricky, it could be fixed when writing documentation. But in order to write documentation it is better to first have test case where this material model is acutally used.

endJunction · 2016-11-08T23:11:20Z

Merge if the config parameter documentation strings are correct.

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

OpenGeoSys development has been moved to GitLab.

See this pull request on GitLab.

endJunction reviewed Nov 3, 2016

View reviewed changes

endJunction approved these changes Nov 3, 2016

View reviewed changes

wenqing force-pushed the perm_Pc_S branch 2 times, most recently from c3c7931 to 4a7c08c Compare November 3, 2016 15:34

endJunction reviewed Nov 3, 2016

View reviewed changes

wenqing force-pushed the perm_Pc_S branch from 0acc833 to e3dff6a Compare November 3, 2016 16:11

wenqing force-pushed the perm_Pc_S branch from e3dff6a to 6cbcbd5 Compare November 7, 2016 07:54

wenqing added the please review label Nov 7, 2016

wenqing added 7 commits November 8, 2016 13:27

[Unsat] Added classes for capillary pressure model

acc88ae

[Unsat] Added Brook-Corey capillary pressure model

2c4ea37

[Unsat] Added van Genuchten capillary pressure model

896933e

[Unsat] Added the creator for capillary pressure models and a mode test

c524fa4

[Unsat] Corrected inline math command of doxygen comment

1b22e4e

[Unsat] Renamed classes and made other minor changes according to Dim…

8f24341

…a's comment

[Unsat] Added a description for a class member and removed an unused …

b492aac

…argument

wenqing force-pushed the perm_Pc_S branch from 6cbcbd5 to b492aac Compare November 8, 2016 12:27

chleh approved these changes Nov 8, 2016

View reviewed changes

wenqing force-pushed the perm_Pc_S branch 2 times, most recently from b7979cf to 76faead Compare November 8, 2016 15:42

[Unsat] Changed the arguments of class constructors and file names, and

688a27a

added model type check in the function body according to the comments by Christoph

wenqing force-pushed the perm_Pc_S branch from 76faead to 688a27a Compare November 8, 2016 15:53

endJunction reviewed Nov 8, 2016

View reviewed changes

endJunction added the please rebase label Nov 8, 2016

endJunction removed the please review label Nov 8, 2016

[Unsat] Corrected keyword documentation.

b4403c8

wenqing merged commit 3b3bfca into ufz:master Nov 9, 2016

wenqing mentioned this pull request Nov 11, 2016

Relative permeability model #1531

Merged

endJunction removed the please rebase label Dec 13, 2016

Capillary pressure models #1517

Capillary pressure models #1517

Conversation

wenqing commented Nov 3, 2016

endJunction Nov 3, 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

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

wenqing commented Nov 3, 2016

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 Nov 3, 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

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 commented Nov 8, 2016

ogsbot commented Jun 19, 2020

endJunction Nov 3, 2016 •

edited

Loading