Merge pull request #1460 from wenqing/porous_medium

Porous medium properties: permeability, storage and porosity 

MaterialLib/CMakeLists.txt

@@ -8,6 +8,10 @@ append_source_files(SOURCES Fluid)
 append_source_files(SOURCES Fluid/Density)
 append_source_files(SOURCES Fluid/Viscosity)
 
+append_source_files(SOURCES PorousMedium/Porosity)
+append_source_files(SOURCES PorousMedium/Storage)
+append_source_files(SOURCES PorousMedium/Permeability)
+
 add_library(MaterialLib ${SOURCES} )
 target_link_libraries(MaterialLib
     BaseLib

MaterialLib/PorousMedium/Permeability/createPermeabilityModel.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file   createPermeabilityModel.cpp
+ *
+ * Created on August 17, 2016, 2:36 PM
+ */
+
+#include "createPermeabilityModel.h"
+
+#include <cassert>
+
+#include "BaseLib/ConfigTree.h"
+
+#include "BaseLib/Error.h"
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+Eigen::MatrixXd createPermeabilityModel(BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{material__porous_medium__permeability__values}
+    auto const values =
+        config.getConfigParameter<std::vector<double>>("values");
+
+    auto const data_size = values.size();
+
+    int dim = -1;
+    switch (data_size)
+    {
+        case 1:
+            dim = 1;
+            break;
+        case 4:
+            dim = 2;
+            break;
+        case 9:
+            dim = 3;
+            break;
+        default:
+        {
+            OGS_FATAL(
+                "Number of values for permeability tensor must be 1, 4 or 9.")
+        }
+    }
+
+    return MathLib::toMatrix(values, dim, dim);
+}
+
+}  // end of namespace
+}  // end of namespace

MaterialLib/PorousMedium/Permeability/createPermeabilityModel.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file   createPermeabilityModel.h
+ *
+ * Created on August 17, 2016, 2:36 PM
+ */
+
+#ifndef OGS_CREATEPERMEABILITYMODEL_H
+#define OGS_CREATEPERMEABILITYMODEL_H
+
+#include <Eigen/Dense>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+/** Create a porosity model
+ *  @param config  ConfigTree object has a tag of <permeability>
+ */
+Eigen::MatrixXd createPermeabilityModel(BaseLib::ConfigTree const& config);
+
+}  // end of namespace
+}  // end of namespace
+
+#endif /* CREATEPERMEABILITYMODEL_H */

MaterialLib/PorousMedium/Porosity/ConstantPorosity.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file   ConstantPorosity.h
+ *
+ * Created on August 16, 2016, 1:03 PM
+ */
+
+#ifndef OGS_CONSTANTPOROSITY_H
+#define OGS_CONSTANTPOROSITY_H
+
+#include "Porosity.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+class ConstantPorosity final : public Porosity
+{
+public:
+    explicit ConstantPorosity(const double value) : _value(value) {}
+
+    /// Get model name.
+    std::string getName() const override { return "Constant porosity"; }
+
+    /**
+     *  Get property value.
+     *  @param variable    A variable with any double type value.
+     *  @param temperature Temperature with any double type value.
+     */
+    double getValue(const double variable, double temperature) const override
+    {
+        (void)variable;
+        (void)temperature;
+        return _value;
+    }
+
+private:
+    const double _value;
+};
+
+}  // end of namespace
+}  // end of namespace
+
+#endif /* CONSTANTPOROSITY_H */

MaterialLib/PorousMedium/Porosity/Porosity.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file:   Porosity.h
+ *
+ * Created on August 16, 2016, 12:53 PM
+ */
+
+#ifndef OGS_POROSITY_H
+#define OGS_POROSITY_H
+
+#include <string>
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+class Porosity
+{
+public:
+    virtual ~Porosity() = default;
+
+    /// Get model name.
+    virtual std::string getName() const = 0;
+
+    /**
+     *  Get property value.
+     *  @param variable    A variable that can be saturation, or an invariant
+     *                     of stress or strain.
+     *  @param temperature Temperature.
+     */
+    virtual double getValue(const double variable,
+                            double temperature) const = 0;
+};
+
+}  // end of namespace
+}  // end of namespace
+
+#endif /* POROSITY_H */

MaterialLib/PorousMedium/Porosity/createPorosityModel.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file:   createPorosityModel.cpp
+ *
+ * Created on August 16, 2016, 1:16 PM
+ */
+
+#include "createPorosityModel.h"
+
+#include "BaseLib/Error.h"
+#include "BaseLib/ConfigTree.h"
+
+#include "Porosity.h"
+#include "ConstantPorosity.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+std::unique_ptr<Porosity> createPorosityModel(BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{material__porous_medium__storage__type}
+    auto const type = config.getConfigParameter<std::string>("type");
+
+    if (type == "Constant")
+        //! \ogs_file_param{material__porous_medium__porosity__Constant_value}
+        return std::unique_ptr<Porosity>(
+            new ConstantPorosity(config.getConfigParameter<double>("value")));
+    else
+    {
+        OGS_FATAL("The porosity type %s is unavailable.\n",
+                  "The available type is \n\tConstant.", type.data());
+    }
+}
+
+}  // end namespace
+}  // end namespace

MaterialLib/PorousMedium/Porosity/createPorosityModel.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file:   createPorosityModel.h
+ *
+ * Created on August 16, 2016, 1:16 PM
+ */
+
+#ifndef OGS_CREATEPOROSITYMODEL_H
+#define OGS_CREATEPOROSITYMODEL_H
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+class Porosity;
+
+/** Create a porosity model
+ *  @param config  ConfigTree object has a tag of <porosity>
+ */
+std::unique_ptr<Porosity> createPorosityModel(
+    BaseLib::ConfigTree const& config);
+
+}  // end namespace
+}  // end namespace
+
+#endif /* CREATEPOROSITYMODEL_H */

MaterialLib/PorousMedium/PorousPropertyHeaders.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * File:   PorousPropertyHeaders.h
+ *
+ * Created on August 18, 2016, 11:12 AM
+ */
+
+#ifndef OGS_POROUSPROPERTYHEADERS_H
+#define OGS_POROUSPROPERTYHEADERS_H
+
+#include "Permeability/createPermeabilityModel.h"
+#include "Porosity/createPorosityModel.h"
+#include "Storage/createStorageModel.h"
+
+#endif /* POROUSPROPERTYHEADERS_H */

MaterialLib/PorousMedium/Storage/ConstantStorage.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file   ConstantStorage.h
+ *
+ * Created on August 16, 2016, 1:03 PM
+ */
+
+#ifndef OGS_CONSTANTSTORAGE_H
+#define OGS_CONSTANTSTORAGE_H
+
+#include "Storage.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+class ConstantStorage final : public Storage
+{
+public:
+    explicit ConstantStorage(const double value) : _value(value) {}
+
+    /// Get model name.
+    std::string getName() const override { return "Constant storage"; }
+
+    /**
+     *  Get storage value.
+     *  @param variable Variable with any double type value
+     */
+    double getValue(const double variable) const override
+    {
+        (void)variable;
+        return _value;
+    }
+
+private:
+    const double _value;
+};
+
+}  // end of namespace
+}  // end of namespace
+
+#endif /* CONSTANTSTORAGE_H */

MaterialLib/PorousMedium/Storage/Storage.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file:   Storage.h
+ *
+ * Created on August 16, 2016, 12:53 PM
+ */
+
+#ifndef OGS_STORAGE_H
+#define OGS_STORAGE_H
+
+#include <string>
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+class Storage
+{
+public:
+    virtual ~Storage() = default;
+    /// Get model name.
+    virtual std::string getName() const = 0;
+
+    /**
+     *  Get property value.
+     *  @param variable A double type variable
+     */
+    virtual double getValue(const double variable) const = 0;
+};
+
+}  // end of namespace
+}  // end of namespace
+
+#endif /* STORAGE_H */

MaterialLib/PorousMedium/Storage/createStorageModel.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file:   createStorageModel.cpp
+ *
+ * Created on August 16, 2016, 1:16 PM
+ */
+
+#include "createStorageModel.h"
+
+#include "BaseLib/Error.h"
+#include "BaseLib/ConfigTree.h"
+
+#include "Storage.h"
+#include "ConstantStorage.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+std::unique_ptr<Storage> createStorageModel(BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{material__porous_medium__storage__type}
+    auto const type = config.getConfigParameter<std::string>("type");
+
+    if (type == "Constant")
+        //! \ogs_file_param{material__porous_medium__storage__constant__value}
+        return std::unique_ptr<Storage>(
+            new ConstantStorage(config.getConfigParameter<double>("value")));
+    else
+    {
+        OGS_FATAL("The storage type %s is unavailable.\n", type.data(),
+                  "The available type is Constant.");
+    }
+}
+
+}  // end namespace
+}  // end namespace

MaterialLib/PorousMedium/Storage/createStorageModel.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file:   createStorageModel.h
+ *
+ * Created on August 16, 2016, 1:16 PM
+ */
+
+#ifndef OGS_CREATESTORAGEMODEL_H
+#define OGS_CREATESTORAGEMODEL_H
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+class Storage;
+/** Create a storage model
+ *  @param config  ConfigTree object has a tag of <storage>
+ */
+std::unique_ptr<Storage> createStorageModel(BaseLib::ConfigTree const& config);
+
+}  // end namespace
+}  // end namespace
+
+#endif /* CREATESTORAGEMODEL_H */

Tests/MaterialLib/TestPermeabilityModel.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ * \file   TestPermeabilityModel.cpp
+ *
+ * Created on August 17, 2016, 4:00 PM
+ */
+#include <gtest/gtest.h>
+
+#include "BaseLib/ConfigTree.h"
+
+#include "TestTools.h"
+
+#include "MaterialLib/PorousMedium/Permeability/createPermeabilityModel.h"
+
+using namespace MaterialLib;
+using namespace MaterialLib::PorousMedium;
+
+using Matrix = Eigen::MatrixXd;
+
+Matrix createTestPermeabilityModel(const char xml[])
+{
+    auto const ptree = readXml(xml);
+    BaseLib::ConfigTree conf(ptree, "", BaseLib::ConfigTree::onerror,
+                             BaseLib::ConfigTree::onwarning);
+    auto const& sub_config = conf.getConfigSubtree("permeability");
+    return MaterialLib::PorousMedium::createPermeabilityModel(sub_config);
+}
+
+class MaterialPermeabilityTest3D : public ::testing::Test
+{
+public:
+    MaterialPermeabilityTest3D() : _dphi(3, 4)
+    {
+        // Mock of the gradient of shape functions
+        _dphi(0, 0) = 1.0;
+        _dphi(0, 1) = 1.0;
+        _dphi(0, 2) = 1.0;
+        _dphi(0, 3) = 1.0;
+        _dphi(1, 0) = 2.0;
+        _dphi(1, 1) = 2.0;
+        _dphi(1, 2) = 2.0;
+        _dphi(1, 3) = 2.0;
+        _dphi(2, 0) = 3.0;
+        _dphi(2, 1) = 3.0;
+        _dphi(2, 2) = 3.0;
+        _dphi(2, 3) = 3.0;
+
+        _trans_dphi = _dphi.transpose();
+    }
+
+protected:
+    Matrix _dphi;
+    Matrix _trans_dphi;
+};
+
+TEST_F(MaterialPermeabilityTest3D, checkAnisotropicPermeability3D)
+{
+    const char xml[] =
+        "<permeability>"
+        "   <values> 2.e-10 0. 0. 0. 3.e-10 0. 0. 0. 4.0e-10 </values> "
+        "</permeability>";
+    Matrix anisK(createTestPermeabilityModel(xml));
+
+    ASSERT_EQ(3, anisK.rows());
+
+    Matrix laplacian = _trans_dphi * anisK * _dphi;
+
+    ASSERT_NEAR(5.e-9, laplacian(0, 0), 1.e-16);
+    ASSERT_NEAR(5.e-9, laplacian(2, 2), 1.e-16);
+}
+
+TEST_F(MaterialPermeabilityTest3D, checkIsotropicPermeability3D)
+{
+    const char xml[] =
+        "<permeability>"
+        "   <values> 2.e-10</values> "
+        "</permeability>";
+    Matrix K(createTestPermeabilityModel(xml));
+
+    ASSERT_EQ(1, K.size());
+
+    const double k = K(0, 0);
+
+    Matrix laplacian = _trans_dphi * k * _dphi;
+
+    ASSERT_NEAR(2.8e-9, laplacian(0, 0), 1.e-16);
+    ASSERT_NEAR(2.8e-9, laplacian(2, 2), 1.e-16);
+}
+
+class MaterialPermeabilityTest2D : public ::testing::Test
+{
+public:
+    MaterialPermeabilityTest2D() : _dphi(2, 3)
+    {
+        // Mock of the gradient of shape functions
+        _dphi(0, 0) = 1.0;
+        _dphi(0, 1) = 1.0;
+        _dphi(0, 2) = 1.0;
+        _dphi(1, 0) = 2.0;
+        _dphi(1, 1) = 2.0;
+        _dphi(1, 2) = 2.0;
+
+        _trans_dphi = _dphi.transpose();
+    }
+
+protected:
+    Matrix _dphi;
+    Matrix _trans_dphi;
+};
+
+TEST_F(MaterialPermeabilityTest2D, checkAnisotropicPermeability2D)
+{
+    const char xml[] =
+        "<permeability>"
+        "   <values> 2.e-10 0. 0. 3.e-10</values> "
+        "</permeability>";
+    Matrix anisK(createTestPermeabilityModel(xml));
+
+    ASSERT_EQ(2, anisK.rows());
+
+    Matrix laplacian = _trans_dphi * anisK * _dphi;
+
+    ASSERT_NEAR(1.4e-9, laplacian(0, 0), 1.e-16);
+    ASSERT_NEAR(1.4e-9, laplacian(2, 2), 1.e-16);
+}
+
+TEST_F(MaterialPermeabilityTest2D, checkIsotropicPermeability2D)
+{
+    const char xml[] =
+        "<permeability>"
+        "   <values> 2.e-10</values> "
+        "</permeability>";
+    Matrix K(createTestPermeabilityModel(xml));
+
+    ASSERT_EQ(1, K.size());
+
+    const double k = K(0, 0);
+
+    Matrix laplacian = _trans_dphi * k * _dphi;
+
+    ASSERT_NEAR(1.0e-9, laplacian(0, 0), 1.e-16);
+    ASSERT_NEAR(1.0e-9, laplacian(2, 2), 1.e-16);
+}

Tests/MaterialLib/TestPorousMediumPorosity.cpp

+/*!
+   \file  TestPorousMediumPorosity.cpp
+   \brief Test the classes for porosity models.
+
+   \copyright
+    Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+               Distributed under a Modified BSD License.
+               See accompanying file LICENSE.txt or
+               http://www.opengeosys.org/project/license
+*/
+
+#include <gtest/gtest.h>
+
+#include <memory>
+
+#include "TestTools.h"
+
+#include "BaseLib/ConfigTree.h"
+
+#include "MaterialLib/PorousMedium/Porosity/ConstantPorosity.h"
+#include "MaterialLib/PorousMedium/Porosity/createPorosityModel.h"
+
+using namespace MaterialLib;
+using namespace MaterialLib::PorousMedium;
+
+std::unique_ptr<Porosity> createTestPorosityModel(const char xml[])
+{
+    auto const ptree = readXml(xml);
+    BaseLib::ConfigTree conf(ptree, "", BaseLib::ConfigTree::onerror,
+                             BaseLib::ConfigTree::onwarning);
+    auto const& sub_config = conf.getConfigSubtree("porosity");
+    return MaterialLib::PorousMedium::createPorosityModel(sub_config);
+}
+
+TEST(Material, checkConstantPorosity)
+{
+    const char xml[] =
+        "<porosity>"
+        "   <type>Constant</type>"
+        "   <value> 0.2 </value> "
+        "</porosity>";
+    auto const n = createTestPorosityModel(xml);
+
+    const double variable = 0.;
+    const double temperature = 0.;
+    ASSERT_EQ(0.2, n->getValue(variable, temperature));
+}

Tests/MaterialLib/TestPorousMediumStorage.cpp

+/*!
+   \file  TestPorousMediumStorage.cpp
+   \brief Test the classes for storage models.
+
+   \copyright
+    Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+               Distributed under a Modified BSD License.
+               See accompanying file LICENSE.txt or
+               http://www.opengeosys.org/project/license
+*/
+
+#include <gtest/gtest.h>
+
+#include <memory>
+
+#include "TestTools.h"
+#include "BaseLib/ConfigTree.h"
+
+#include "MaterialLib/PorousMedium/Storage/ConstantStorage.h"
+#include "MaterialLib/PorousMedium/Storage/createStorageModel.h"
+
+using namespace MaterialLib;
+using namespace MaterialLib::PorousMedium;
+
+std::unique_ptr<Storage> createTestStorageModel(const char xml[])
+{
+    auto const ptree = readXml(xml);
+    BaseLib::ConfigTree conf(ptree, "", BaseLib::ConfigTree::onerror,
+                             BaseLib::ConfigTree::onwarning);
+    auto const& sub_config = conf.getConfigSubtree("storage");
+    return MaterialLib::PorousMedium::createStorageModel(sub_config);
+}
+
+TEST(Material, checkConstantStorage)
+{
+    const char xml[] =
+        "<storage>"
+        "   <type>Constant</type>"
+        "   <value> 1.e-4 </value> "
+        "</storage>";
+    auto const eta = createTestStorageModel(xml);
+
+    const double var = 0;
+    ASSERT_EQ(1.e-4, eta->getValue(var));
+}
+