diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index a035a3fada6418ef29c00aaed21f36225165553d..19c2160f56fd45cdf02ed3622b3bb7953171d71f 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -35,6 +35,7 @@
#include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h"
#include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h"
#include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
+#include "ProcessLib/LiquidFlow/CreateLiquidFlowProcess.h"
#include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
#include "ProcessLib/SmallDeformationWithLIE/CreateSmallDeformationProcess.h"
#include "ProcessLib/TES/CreateTESProcess.h"
@@ -292,6 +293,12 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
_process_variables, _parameters, integration_order,
process_config, project_directory, output_directory);
}
+ else if (type == "LIQUID_FLOW")
+ {
+ process = ProcessLib::LiquidFlow::createLiquidFlowProcess(
+ *_mesh_vec[0], std::move(jacobian_assembler), _process_variables,
+ _parameters, integration_order, process_config);
+ }
else if (type == "TES")
{
process = ProcessLib::TES::createTESProcess(
diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index 4c4fe0f0a1ade1be5403baf6c82d36b61f54287f..174b023608b0858e1e404cd8b462b6c3a64fbf60 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -574,6 +574,52 @@ if(NOT OGS_USE_MPI)
single_joint_inside_expected_pcs_0_ts_1_t_1.000000.vtu single_joint_inside_pcs_0_ts_1_t_1.000000.vtu displacement_jump1 displacement_jump1
)
+ # Liquid flow
+ AddTest(
+ NAME LiquidFlow_LineDirichletNeumannBC
+ PATH Parabolic/LiquidFlow/LineDirichletNeumannBC
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS line_dirichlet_neumannBC.prj
+ WRAPPER time
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-8
+ DIFF_DATA
+ sat1D.vtu sat_1D_pcs_0_ts_1_t_1.000000.vtu AnalyticSolution pressure
+ )
+ AddTest(
+ NAME LiquidFlow_PressureBCatCornerOfAnisotropicSquare
+ PATH Parabolic/LiquidFlow/PressureBCatCornerOfAnisotropicSquare
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS pressureBC_at_corner_of_anisotropic_square.prj
+ WRAPPER time
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-8
+ DIFF_DATA
+ mesh2D.vtu sat_2D_lflow_pcs_0_ts_1_t_1.000000.vtu OGS5_Results pressure
+ )
+ AddTest(
+ NAME LiquidFlow_GravityDriven
+ PATH Parabolic/LiquidFlow/GravityDriven
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS gravity_driven.prj
+ WRAPPER time
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-8
+ DIFF_DATA
+ mesh2D.vtu gravity_driven_pcs_0_ts_1_t_1.000000.vtu AnalyticSolution pressure
+ )
+ AddTest(
+ NAME LiquidFlow_AxisymTheis
+ PATH Parabolic/LiquidFlow/AxiSymTheis
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS axisym_theis.prj
+ WRAPPER time
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-8
+ DIFF_DATA
+ axisym_theis.vtu liquid_pcs_pcs_0_ts_30_t_1728.000000.vtu OGS5_pressure pressure
+ )
+
else()
# MPI groundwater flow tests
AddTest(
@@ -782,4 +828,51 @@ else()
tes_zeolite_discharge_large_ts_28_t_1_000000.vtu tes_zeolite_discharge_large_pcs_0_ts_28_t_1_000000_0.vtu v_mass_frac v_mass_frac
# tes_zeolite_discharge_large_ts_28_t_1_0.vtu solid_density solid_density
)
+
+ # Liquid flow
+ AddTest(
+ NAME LiquidFlow_LineDirichletNeumannBC
+ PATH Parabolic/LiquidFlow/LineDirichletNeumannBC
+ EXECUTABLE_ARGS line_dirichlet_neumannBC.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-8
+ DIFF_DATA
+ sat1D.vtu sat_1D_pcs_0_ts_1_t_1.000000.vtu AnalyticSolution pressure
+ )
+ AddTest(
+ NAME LiquidFlow_GravityDriven
+ PATH Parabolic/LiquidFlow/GravityDriven
+ EXECUTABLE_ARGS gravity_driven.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-8
+ DIFF_DATA
+ mesh2D.vtu gravity_driven_pcs_0_ts_1_t_1.000000.vtu AnalyticSolution pressure
+ )
+ AddTest(
+ NAME LiquidFlow_PressureBCatCornerOfAnisotropicSquare
+ PATH Parabolic/LiquidFlow/PressureBCatCornerOfAnisotropicSquare
+ EXECUTABLE_ARGS pressureBC_at_corner_of_anisotropic_square.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-8
+ DIFF_DATA
+ mesh2D.vtu sat_2D_lflow_pcs_0_ts_1_t_1.000000.vtu OGS5_Results pressure
+ )
+ AddTest(
+ NAME LiquidFlow_AxisymTheis
+ PATH Parabolic/LiquidFlow/AxiSymTheis
+ EXECUTABLE_ARGS axisym_theis.prj
+ WRAPPER mpirun
+ WRAPPER_ARGS -np 1
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-8
+ DIFF_DATA
+ axisym_theis.vtu liquid_pcs_pcs_0_ts_30_t_1728.000000.vtu OGS5_pressure pressure
+ )
+
endif()
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index 7be84b171a279277277b65758449704dc51ea699..0dd7a458c2bf1ea8f58e376dfebf9662a5c3cf53 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -14,6 +14,8 @@ APPEND_SOURCE_FILES(SOURCES Parameter)
add_subdirectory(GroundwaterFlow)
APPEND_SOURCE_FILES(SOURCES GroundwaterFlow)
+APPEND_SOURCE_FILES(SOURCES LiquidFlow)
+
APPEND_SOURCE_FILES(SOURCES SmallDeformation)
APPEND_SOURCE_FILES(SOURCES SmallDeformationWithLIE)
diff --git a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..78db4b6026b1f452c3a9d8813e20cd6c83e1194a
--- /dev/null
+++ b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.cpp
@@ -0,0 +1,105 @@
+/**
+ * \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 CreateLiquidFlowProcess.cpp
+ *
+ * Created on August 19, 2016, 1:30 PM
+ */
+#include "CreateLiquidFlowProcess.h"
+
+#include <algorithm>
+
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+#include "LiquidFlowProcess.h"
+
+namespace ProcessLib
+{
+namespace LiquidFlow
+{
+std::unique_ptr<Process> createLiquidFlowProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{process__type}
+ config.checkConfigParameter("type", "LIQUID_FLOW");
+
+ DBUG("Create LiquidFlowProcess.");
+
+ // Process variable.
+ auto process_variables = findProcessVariables(
+ variables, config,
+ {//! \ogs_file_param_special{process__LIQUID_FLOW__process_variables__process_variable}
+ "process_variable"});
+
+ SecondaryVariableCollection secondary_variables;
+
+ NumLib::NamedFunctionCaller named_function_caller({"LiquidFlow_pressure"});
+
+ ProcessLib::parseSecondaryVariables(config, secondary_variables,
+ named_function_caller);
+
+ // Get the gravity vector for the Darcy velocity
+ //! \ogs_file_param{process__LIQUID_FLOW__darcy_gravity}
+ auto const& darcy_g_config = config.getConfigSubtree("darcy_gravity");
+ const int gravity_axis_id_input =
+ //! \ogs_file_param_special{process__LIQUID_FLOW__darcy_gravity_axis_id}
+ darcy_g_config.getConfigParameter<int>("axis_id");
+ assert(gravity_axis_id_input < static_cast<int>(mesh.getDimension()));
+ const double g =
+ //! \ogs_file_param_special{process__LIQUID_FLOW__darcy_gravity_g}
+ darcy_g_config.getConfigParameter<double>("g");
+ assert(g >= 0.);
+ const int gravity_axis_id = (g == 0.) ? -1 : gravity_axis_id_input;
+
+ //! \ogs_file_param{process__LIQUID_FLOW__material_property}
+ auto const& mat_config = config.getConfigSubtree("material_property");
+
+ auto const& mat_ids =
+ mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+ if (mat_ids)
+ {
+ INFO("The liquid flow is in heterogeneous porous media.");
+ const bool has_material_ids = true;
+ return std::unique_ptr<Process>{new LiquidFlowProcess{
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(secondary_variables),
+ std::move(named_function_caller), *mat_ids, has_material_ids,
+ gravity_axis_id, g, mat_config}};
+ }
+ else
+ {
+ INFO("The liquid flow is in homogeneous porous media.");
+
+ MeshLib::Properties dummy_property;
+ // For a reference argument of LiquidFlowProcess(...).
+ auto const& dummy_property_vector =
+ dummy_property.createNewPropertyVector<int>(
+ "MaterialIDs", MeshLib::MeshItemType::Cell, 1);
+
+ // Since dummy_property_vector is only visible in this function,
+ // the following constant, has_material_ids, is employed to indicate
+ // that material_ids does not exist.
+ const bool has_material_ids = false;
+
+ return std::unique_ptr<Process>{new LiquidFlowProcess{
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(secondary_variables),
+ std::move(named_function_caller), *dummy_property_vector,
+ has_material_ids, gravity_axis_id, g, mat_config}};
+ }
+}
+
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.h b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..6e199f53db994f5611d7155901847e2c96f25e2c
--- /dev/null
+++ b/ProcessLib/LiquidFlow/CreateLiquidFlowProcess.h
@@ -0,0 +1,33 @@
+/**
+ * \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 CreateLiquidFlowProcess.h
+ *
+ * Created on August 19, 2016, 1:30 PM
+ */
+
+#ifndef OGS_CREATELIQUIDFLOWPROCESS_H
+#define OGS_CREATELIQUIDFLOWPROCESS_H
+
+#include <memory>
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace LiquidFlow
+{
+std::unique_ptr<Process> createLiquidFlowProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<ProcessVariable> const& variables,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ BaseLib::ConfigTree const& config);
+} // end of namespace
+} // end of namespace
+
+#endif /* CREATELIQUIDFLOWPROCESS_H */
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..34c42024fdb6badc912870991fdd1e4e6221b3b1
--- /dev/null
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler-impl.h
@@ -0,0 +1,156 @@
+/**
+ * \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 LiquidFlowLocalAssembler.cpp
+ *
+ * Created on August 19, 2016, 2:28 PM
+ */
+
+#ifndef OGS_LIQUIDFLOWLOCALASSEMBLER_IMPL_H
+#define OGS_LIQUIDFLOWLOCALASSEMBLER_IMPL_H
+
+#include "LiquidFlowLocalAssembler.h"
+
+#include "NumLib/Function/Interpolation.h"
+
+namespace ProcessLib
+{
+namespace LiquidFlow
+{
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ assemble(double const t, std::vector<double> const& local_x,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data)
+{
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+ _material_properties.setMaterialID(pos);
+
+ const Eigen::MatrixXd& perm =
+ _material_properties.getPermeability(t, pos, _element.getDimension());
+ // Note: For Inclined 1D in 2D/3D or 2D element in 3D, the first item in
+ // the assert must be changed to perm.rows() == _element->getDimension()
+ assert(perm.rows() == GlobalDim || perm.rows() == 1);
+
+ if (perm.size() == 1) // isotropic or 1D problem.
+ local_assemble<IsotropicLaplacianAndGravityTermCalculator>(
+ t, local_x, local_M_data, local_K_data, local_b_data, pos, perm);
+ else
+ local_assemble<AnisotropicLaplacianAndGravityTermCalculator>(
+ t, local_x, local_M_data, local_K_data, local_b_data, pos, perm);
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+template <typename LaplacianAndGravityTermCalculator>
+void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ local_assemble(double const t, std::vector<double> const& local_x,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data,
+ SpatialPosition const& pos, Eigen::MatrixXd const& perm)
+{
+ auto const local_matrix_size = local_x.size();
+ assert(local_matrix_size == ShapeFunction::NPOINTS);
+
+ auto local_M = MathLib::createZeroedMatrix<NodalMatrixType>(
+ local_M_data, local_matrix_size, local_matrix_size);
+ auto local_K = MathLib::createZeroedMatrix<NodalMatrixType>(
+ local_K_data, local_matrix_size, local_matrix_size);
+ auto local_b = MathLib::createZeroedVector<NodalVectorType>(
+ local_b_data, local_matrix_size);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ // TODO: The following two variables should be calculated inside the
+ // the integration loop for non-constant porosity and storage models.
+ double porosity_variable = 0.;
+ double storage_variable = 0.;
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& sm = _shape_matrices[ip];
+ auto const& wp = _integration_method.getWeightedPoint(ip);
+
+ double p = 0.;
+ NumLib::shapeFunctionInterpolate(local_x, sm.N, p);
+ // TODO : compute _temperature from the heat transport pcs
+
+ const double integration_factor =
+ sm.integralMeasure * sm.detJ * wp.getWeight();
+
+ // Assemble mass matrix, M
+ local_M.noalias() +=
+ _material_properties.getMassCoefficient(
+ t, pos, porosity_variable, storage_variable, p, _temperature) *
+ sm.N.transpose() * sm.N * integration_factor;
+
+ // Compute density:
+ const double rho_g =
+ _material_properties.getLiquidDensity(p, _temperature) *
+ _gravitational_acceleration;
+ // Compute viscosity:
+ const double mu = _material_properties.getViscosity(p, _temperature);
+
+ // Assemble Laplacian, K, and RHS by the gravitational term
+ LaplacianAndGravityTermCalculator::calculate(
+ local_K, local_b, sm, perm, integration_factor, mu, rho_g,
+ _gravitational_axis_id);
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ IsotropicLaplacianAndGravityTermCalculator::calculate(
+ Eigen::Map<NodalMatrixType>& local_K,
+ Eigen::Map<NodalVectorType>& local_b, ShapeMatrices const& sm,
+ Eigen::MatrixXd const& perm, double const integration_factor,
+ double const mu, double const rho_g, int const gravitational_axis_id)
+{
+ const double K = perm(0, 0) / mu;
+ const double fac = K * integration_factor;
+ local_K.noalias() += fac * sm.dNdx.transpose() * sm.dNdx;
+
+ if (gravitational_axis_id >= 0)
+ {
+ // Note: Since perm, K, is a scalar number in this function,
+ // (gradN)*K*V becomes K*(grad N)*V. With the gravity vector of V,
+ // the simplification of (grad N)*V is the gravitational_axis_id th
+ // column of the transpose of (grad N) multiplied with rho_g.
+ local_b.noalias() -=
+ fac * sm.dNdx.transpose().col(gravitational_axis_id) * rho_g;
+ }
+}
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void LiquidFlowLocalAssembler<ShapeFunction, IntegrationMethod, GlobalDim>::
+ AnisotropicLaplacianAndGravityTermCalculator::calculate(
+ Eigen::Map<NodalMatrixType>& local_K,
+ Eigen::Map<NodalVectorType>& local_b, ShapeMatrices const& sm,
+ Eigen::MatrixXd const& perm, double const integration_factor,
+ double const mu, double const rho_g, int const gravitational_axis_id)
+{
+ const double fac = integration_factor / mu;
+ local_K.noalias() += fac * sm.dNdx.transpose() * perm * sm.dNdx;
+ if (gravitational_axis_id >= 0)
+ {
+ // Note: perm * gravity_vector = rho_g * perm.col(gravitational_axis_id)
+ // i.e the result is the gravitational_axis_id th column of
+ // perm multiplied with rho_g.
+ local_b.noalias() -=
+ fac * rho_g * sm.dNdx.transpose() * perm.col(gravitational_axis_id);
+ }
+}
+} // end of namespace
+} // end of namespace
+
+#endif /* OGS_LIQUIDFLOWLOCALASSEMBLER_IMPL_H */
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..6b34deaa09d3a3d943147f8d9bb0672be49daf96
--- /dev/null
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -0,0 +1,178 @@
+/**
+ * \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 LiquidFlowLocalAssembler.h
+ *
+ * Created on August 19, 2016, 2:28 PM
+ */
+
+#ifndef OGS_LIQUIDFLOWLOCALASSEMBLER_H
+#define OGS_LIQUIDFLOWLOCALASSEMBLER_H
+
+#include <vector>
+
+#include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/Extrapolation/ExtrapolatableElement.h"
+#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
+#include "ProcessLib/LocalAssemblerInterface.h"
+#include "ProcessLib/LocalAssemblerTraits.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Utils/InitShapeMatrices.h"
+
+#include "LiquidFlowMaterialProperties.h"
+
+namespace ProcessLib
+{
+namespace LiquidFlow
+{
+const unsigned NUM_NODAL_DOF = 1;
+
+class LiquidFlowLocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+public:
+ virtual std::vector<double> const& getIntPtDarcyVelocityX(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityY(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtDarcyVelocityZ(
+ std::vector<double>& /*cache*/) const = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class LiquidFlowLocalAssembler : public LiquidFlowLocalAssemblerInterface
+{
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+ using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits<
+ ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
+
+ using NodalMatrixType = typename LocalAssemblerTraits::LocalMatrix;
+ using NodalVectorType = typename LocalAssemblerTraits::LocalVector;
+ using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+
+public:
+ LiquidFlowLocalAssembler(MeshLib::Element const& element,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ int const gravitational_axis_id,
+ double const gravitational_acceleration,
+ LiquidFlowMaterialProperties& material_propertries)
+ : _element(element),
+ _integration_method(integration_order),
+ _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ element, is_axially_symmetric, _integration_method)),
+ _gravitational_axis_id(gravitational_axis_id),
+ _gravitational_acceleration(gravitational_acceleration),
+ _material_properties(material_propertries)
+ {
+ }
+
+ void assemble(double const t, std::vector<double> const& local_x,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data) override;
+
+ Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
+ const unsigned integration_point) const override
+ {
+ auto const& N = _shape_matrices[integration_point].N;
+
+ // assumes N is stored contiguously in memory
+ return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocityX(
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 0);
+ return _darcy_velocities[0];
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocityY(
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 1);
+ return _darcy_velocities[1];
+ }
+
+ std::vector<double> const& getIntPtDarcyVelocityZ(
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_darcy_velocities.size() > 2);
+ return _darcy_velocities[2];
+ }
+
+private:
+ MeshLib::Element const& _element;
+
+ IntegrationMethod const _integration_method;
+ std::vector<ShapeMatrices> _shape_matrices;
+
+ std::vector<std::vector<double>> _darcy_velocities =
+ std::vector<std::vector<double>>(
+ GlobalDim,
+ std::vector<double>(_integration_method.getNumberOfPoints()));
+
+ /**
+ * Calculator of the Laplacian and the gravity term for anisotropic
+ * permeability tensor
+ */
+ struct AnisotropicLaplacianAndGravityTermCalculator
+ {
+ static void calculate(Eigen::Map<NodalMatrixType>& local_K,
+ Eigen::Map<NodalVectorType>& local_b,
+ ShapeMatrices const& sm,
+ Eigen::MatrixXd const& perm,
+ double const integration_factor, double const mu,
+ double const rho_g,
+ int const gravitational_axis_id);
+ };
+
+ /**
+ * Calculator of the Laplacian and the gravity term for isotropic
+ * permeability tensor
+ */
+ struct IsotropicLaplacianAndGravityTermCalculator
+ {
+ static void calculate(Eigen::Map<NodalMatrixType>& local_K,
+ Eigen::Map<NodalVectorType>& local_b,
+ ShapeMatrices const& sm,
+ Eigen::MatrixXd const& perm,
+ double const integration_factor, double const mu,
+ double const rho_g,
+ int const gravitational_axis_id);
+ };
+
+ template <typename LaplacianAndGravityTermCalculator>
+ void local_assemble(double const t, std::vector<double> const& local_x,
+ std::vector<double>& local_M_data,
+ std::vector<double>& local_K_data,
+ std::vector<double>& local_b_data,
+ SpatialPosition const& pos,
+ Eigen::MatrixXd const& perm);
+
+ const int _gravitational_axis_id;
+ const double _gravitational_acceleration;
+ LiquidFlowMaterialProperties& _material_properties;
+ double _temperature;
+};
+
+} // end of namespace
+} // end of namespace
+
+#include "LiquidFlowLocalAssembler-impl.h"
+
+#endif /* LIQUIDFLOWLOCALASSEMBLER_H */
diff --git a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..063a1cfb6cd8a11f08e883c46db9d9c092f43788
--- /dev/null
+++ b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.cpp
@@ -0,0 +1,140 @@
+/**
+ * \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 LiquidFlowMaterialProperties.cpp
+ *
+ * Created on August 18, 2016, 11:49 AM
+ */
+
+#include "LiquidFlowMaterialProperties.h"
+
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/reorderVector.h"
+
+#include "MeshLib/PropertyVector.h"
+
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Parameter/SpatialPosition.h"
+
+#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+
+namespace ProcessLib
+{
+namespace LiquidFlow
+{
+LiquidFlowMaterialProperties::LiquidFlowMaterialProperties(
+ BaseLib::ConfigTree const& config,
+ bool const has_material_ids,
+ MeshLib::PropertyVector<int> const& material_ids)
+ : _has_material_ids(has_material_ids), _material_ids(material_ids)
+{
+ DBUG("Reading material properties of liquid flow process.");
+
+ //! \ogs_file_param{prj__material_property__fluid}
+ auto const& fluid_config = config.getConfigSubtree("fluid");
+
+ // Get fluid properties
+ //! \ogs_file_param{prj__material_property__fluid__density}
+ auto const& rho_conf = fluid_config.getConfigSubtree("density");
+ _liquid_density = MaterialLib::Fluid::createFluidDensityModel(rho_conf);
+ //! \ogs_file_param{prj__material_property__fluid__viscosity}
+ auto const& mu_conf = fluid_config.getConfigSubtree("viscosity");
+ _viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf);
+
+ // Get porous properties
+ std::vector<int> mat_ids;
+ //! \ogs_file_param{prj__material_property__porous_medium}
+ auto const& poro_config = config.getConfigSubtree("porous_medium");
+ //! \ogs_file_param{prj__material_property__porous_medium__porous_medium}
+ for (auto const& conf : poro_config.getConfigSubtreeList("porous_medium"))
+ {
+ //! \ogs_file_attr{prj__material_property__porous_medium__porous_medium__id}
+ auto const id = conf.getConfigAttributeOptional<int>("id");
+ mat_ids.push_back(*id);
+
+ //! \ogs_file_param{prj__material_property__porous_medium__porous_medium__permeability}
+ auto const& perm_conf = conf.getConfigSubtree("permeability");
+ _intrinsic_permeability_models.emplace_back(
+ MaterialLib::PorousMedium::createPermeabilityModel(perm_conf));
+
+ //! \ogs_file_param{prj__material_property__porous_medium__porous_medium__porosity}
+ auto const& poro_conf = conf.getConfigSubtree("porosity");
+ auto n = MaterialLib::PorousMedium::createPorosityModel(poro_conf);
+ _porosity_models.emplace_back(std::move(n));
+
+ //! \ogs_file_param{prj__material_property__porous_medium__porous_medium__storage}
+ auto const& stora_conf = conf.getConfigSubtree("storage");
+ auto beta = MaterialLib::PorousMedium::createStorageModel(stora_conf);
+ _storage_models.emplace_back(std::move(beta));
+ }
+
+ BaseLib::reorderVector(_intrinsic_permeability_models, mat_ids);
+ BaseLib::reorderVector(_porosity_models, mat_ids);
+ BaseLib::reorderVector(_storage_models, mat_ids);
+}
+
+void LiquidFlowMaterialProperties::setMaterialID(const SpatialPosition& pos)
+{
+ if (!_has_material_ids)
+ {
+ _current_material_id = 0;
+ return;
+ }
+
+ assert(pos.getElementID().get() < _material_ids.size());
+ _current_material_id = _material_ids[pos.getElementID().get()];
+}
+
+double LiquidFlowMaterialProperties::getLiquidDensity(const double p,
+ const double T) const
+{
+ ArrayType vars;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::pl)] = p;
+ return _liquid_density->getValue(vars);
+}
+
+double LiquidFlowMaterialProperties::getViscosity(const double p,
+ const double T) const
+{
+ ArrayType vars;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::pl)] = p;
+ return _viscosity->getValue(vars);
+}
+
+double LiquidFlowMaterialProperties::getMassCoefficient(
+ const double /*t*/, const SpatialPosition& /*pos*/, const double p,
+ const double T, const double porosity_variable,
+ const double storage_variable) const
+{
+ ArrayType vars;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::T)] = T;
+ vars[static_cast<int>(MaterialLib::Fluid::PropertyVariableType::pl)] = p;
+ const double drho_dp = _liquid_density->getdValue(
+ vars, MaterialLib::Fluid::PropertyVariableType::pl);
+ const double rho = _liquid_density->getValue(vars);
+ assert(rho > 0.);
+
+ const double porosity =
+ _porosity_models[_current_material_id]->getValue(porosity_variable, T);
+ const double storage =
+ _storage_models[_current_material_id]->getValue(storage_variable);
+ return porosity * drho_dp / rho + storage;
+}
+
+Eigen::MatrixXd const& LiquidFlowMaterialProperties::getPermeability(
+ const double /*t*/, const SpatialPosition& /*pos*/, const int /*dim*/) const
+{
+ return _intrinsic_permeability_models[_current_material_id];
+}
+
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..bf26617da8b62430396231f363247656d3f70291
--- /dev/null
+++ b/ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h
@@ -0,0 +1,109 @@
+/**
+ * \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 LiquidFlowMaterialProperties.h
+ *
+ * Created on August 18, 2016, 11:03 AM
+ */
+
+#ifndef OGS_LIQUIDFLOWMATERIALPROPERTIES_H
+#define OGS_LIQUIDFLOWMATERIALPROPERTIES_H
+
+#include <memory>
+
+#include "MaterialLib/Fluid/FluidPropertyHeaders.h"
+#include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+class Porosity;
+class Storage;
+}
+}
+
+namespace MeshLib
+{
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
+namespace ProcessLib
+{
+class SpatialPosition;
+
+namespace LiquidFlow
+{
+class LiquidFlowMaterialProperties
+{
+public:
+ typedef MaterialLib::Fluid::FluidProperty::ArrayType ArrayType;
+
+ LiquidFlowMaterialProperties(
+ BaseLib::ConfigTree const& config,
+ bool const has_material_ids,
+ MeshLib::PropertyVector<int> const& material_ids);
+
+ void setMaterialID(const SpatialPosition& pos);
+
+ /**
+ * \brief Compute the coefficient of the mass term by
+ * \f[
+ * n \frac{partial \rho_l}{\partial p} + \beta_s
+ * \f]
+ * where \f$n\f$ is the porosity, \f$rho_l\f$ is the liquid density,
+ * \f$bata_s\f$ is the storage.
+ * \param t Time.
+ * \param pos Position of element.
+ * \param p Pressure value.
+ * \param T Temperature value.
+ * \param porosity_variable The first variable for porosity model, and it
+ * passes a double type value that could be
+ * saturation, and invariant of stress or strain.
+ * \param storage_variable Variable for storage model.
+ */
+ double getMassCoefficient(const double t, const SpatialPosition& pos,
+ const double p, const double T,
+ const double porosity_variable,
+ const double storage_variable) const;
+
+ Eigen::MatrixXd const& getPermeability(const double t,
+ const SpatialPosition& pos,
+ const int dim) const;
+
+ double getLiquidDensity(const double p, const double T) const;
+
+ double getViscosity(const double p, const double T) const;
+
+private:
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_density;
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty> _viscosity;
+
+ /// A flag to indicate whether the reference member, _material_ids,
+ /// is not assigned.
+ const bool _has_material_ids;
+ /** Use porous medium models for different material zones.
+ * Material IDs must be given as mesh element properties.
+ */
+ MeshLib::PropertyVector<int> const& _material_ids;
+
+ int _current_material_id = 0;
+ std::vector<Eigen::MatrixXd> _intrinsic_permeability_models;
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Porosity>>
+ _porosity_models;
+ std::vector<std::unique_ptr<MaterialLib::PorousMedium::Storage>>
+ _storage_models;
+
+ // Note: For the statistical data of porous media, they could be read from
+ // vtu files directly. This can be done by using property vectors directly.
+ // Such property vectors will be added here if they are needed.
+};
+
+} // end of namespace
+} // end of namespace
+#endif /* LIQUIDFLOWMATERIALPROPERTIES_H */
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..603d30d776f173bae0aeacfd59fc543cb988d7a6
--- /dev/null
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.cpp
@@ -0,0 +1,115 @@
+/**
+ * \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 LiquidFlowProcess.cpp
+ *
+ * Created on August 19, 2016, 1:38 PM
+ */
+
+#include "LiquidFlowProcess.h"
+
+#include <cassert>
+
+#include "MeshLib/PropertyVector.h"
+
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "LiquidFlowLocalAssembler.h"
+
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+
+namespace ProcessLib
+{
+namespace LiquidFlow
+{
+LiquidFlowProcess::LiquidFlowProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ std::vector<std::reference_wrapper<ProcessVariable>>&& process_variables,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller,
+ MeshLib::PropertyVector<int> const& material_ids,
+ bool const has_material_ids,
+ int const gravitational_axis_id,
+ double const gravitational_acceleration,
+ BaseLib::ConfigTree const& config)
+ : Process(mesh, std::move(jacobian_assembler), parameters,
+ integration_order, std::move(process_variables),
+ std::move(secondary_variables), std::move(named_function_caller)),
+ _gravitational_axis_id(gravitational_axis_id),
+ _gravitational_acceleration(gravitational_acceleration),
+ _material_properties(
+ LiquidFlowMaterialProperties(config, has_material_ids, material_ids))
+{
+ DBUG("Create Liquid flow process.");
+}
+
+void LiquidFlowProcess::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ ProcessLib::createLocalAssemblers<LiquidFlowLocalAssembler>(
+ mesh.getDimension(), mesh.getElements(), dof_table, _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, _gravitational_axis_id,
+ _gravitational_acceleration, _material_properties);
+
+ _secondary_variables.addSecondaryVariable(
+ "darcy_velocity_x", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityX));
+
+ if (mesh.getDimension() > 1)
+ {
+ _secondary_variables.addSecondaryVariable(
+ "darcy_velocity_y", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityY));
+ }
+ if (mesh.getDimension() > 2)
+ {
+ _secondary_variables.addSecondaryVariable(
+ "darcy_velocity_z", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &LiquidFlowLocalAssemblerInterface::getIntPtDarcyVelocityZ));
+ }
+}
+
+void LiquidFlowProcess::assembleConcreteProcess(const double t,
+ GlobalVector const& x,
+ GlobalMatrix& M,
+ GlobalMatrix& K,
+ GlobalVector& b)
+{
+ DBUG("Assemble LiquidFlowProcess.");
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assemble, _local_assemblers,
+ *_local_to_global_index_map, t, x, M, K, b);
+}
+
+void LiquidFlowProcess::assembleWithJacobianConcreteProcess(
+ const double t, GlobalVector const& x, GlobalVector const& xdot,
+ const double dxdot_dx, const double dx_dx, GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b, GlobalMatrix& Jac)
+{
+ DBUG("AssembleWithJacobian LiquidFlowProcess.");
+
+ // Call global assembler for each local assembly item.
+ GlobalExecutor::executeMemberDereferenced(
+ _global_assembler, &VectorMatrixAssembler::assembleWithJacobian,
+ _local_assemblers, *_local_to_global_index_map, t, x, xdot, dxdot_dx,
+ dx_dx, M, K, b, Jac);
+}
+
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/LiquidFlow/LiquidFlowProcess.h b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..ff1a5aac2052550c6cfd6c6cf592925b79c09ffc
--- /dev/null
+++ b/ProcessLib/LiquidFlow/LiquidFlowProcess.h
@@ -0,0 +1,98 @@
+/**
+ * \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 LiquidFlowProcess.h
+ *
+ * Created on August 19, 2016, 1:38 PM
+ */
+
+#ifndef OGS_LIQUIDFLOWPROCESS_H
+#define OGS_LIQUIDFLOWPROCESS_H
+
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/Process.h"
+
+#include "LiquidFlowMaterialProperties.h"
+#include "LiquidFlowLocalAssembler.h"
+
+namespace MeshLib
+{
+class Element;
+class Mesh;
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
+namespace ProcessLib
+{
+namespace LiquidFlow
+{
+/**
+ * \brief A class to simulate the liquid flow process in porous media described
+ * by
+ *
+ * \f[
+ * \frac{\partial n \rho_l}{\partial T} \frac{\partial T}{\partial t}/\rho_l
+ * + (\frac{\partial n \rho_l}{\partial p}/\rho_l + \beta_s)
+ * \frac{\partial p}{\partial t}
+ * -\nabla (\frac{K}{\mu}(\nabla p + \rho_l g \nabla z) ) = Q
+ * \f]
+ * where
+ * \f{eqnarray*}{
+ * &p:& \mbox{pore pressure,}\\
+ * &T: & \mbox{Temperature,}\\
+ * &\rho_l:& \mbox{liquid density,}\\
+ * &\beta_s:& \mbox{specific storage,}\\
+ * &K:& \mbox{permeability,}\\
+ * &\mu:& \mbox{viscosity,}\\
+ * \f}
+ */
+class LiquidFlowProcess final : public Process
+{
+public:
+ LiquidFlowProcess(
+ MeshLib::Mesh& mesh,
+ std::unique_ptr<AbstractJacobianAssembler>&& jacobian_assembler,
+ std::vector<std::unique_ptr<ParameterBase>> const& parameters,
+ unsigned const integration_order,
+ std::vector<std::reference_wrapper<ProcessVariable>>&&
+ process_variables,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller,
+ MeshLib::PropertyVector<int> const& material_ids,
+ bool const has_material_ids,
+ int const gravitational_axis_id,
+ double const gravitational_acceleration,
+ BaseLib::ConfigTree const& config);
+
+ bool isLinear() const override { return true; }
+private:
+ void initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh, unsigned const integration_order) override;
+
+ void assembleConcreteProcess(const double t, GlobalVector const& x,
+ GlobalMatrix& M, GlobalMatrix& K,
+ GlobalVector& b) override;
+
+ void assembleWithJacobianConcreteProcess(
+ const double t, GlobalVector const& x, GlobalVector const& xdot,
+ const double dxdot_dx, const double dx_dx, GlobalMatrix& M,
+ GlobalMatrix& K, GlobalVector& b, GlobalMatrix& Jac) override;
+
+ const int _gravitational_axis_id;
+ const double _gravitational_acceleration;
+ LiquidFlowMaterialProperties _material_properties;
+
+ std::vector<std::unique_ptr<LiquidFlowLocalAssemblerInterface>>
+ _local_assemblers;
+};
+
+} // end of namespace
+} // end of namespace
+
+#endif /* LIQUIDFLOWPROCESS_H */
diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt
index 18735e31f903a258203f56e7611765d5415df94b..cb082141837322f0e774d2e63dbf3dea16ad2cd1 100644
--- a/Tests/CMakeLists.txt
+++ b/Tests/CMakeLists.txt
@@ -18,6 +18,7 @@ APPEND_SOURCE_FILES(TEST_SOURCES MeshLib)
APPEND_SOURCE_FILES(TEST_SOURCES MeshGeoToolsLib)
APPEND_SOURCE_FILES(TEST_SOURCES NumLib)
APPEND_SOURCE_FILES(TEST_SOURCES ProcessLib)
+APPEND_SOURCE_FILES(TEST_SOURCES ProcessLib/LiquidFlow)
if(QT4_FOUND)
APPEND_SOURCE_FILES(TEST_SOURCES FileIO_Qt)
diff --git a/Tests/Data b/Tests/Data
index 6de9e6af6d103d65054b2f2b16f6d6b49a53aeaf..d6ca72279bdc1275b7ccd98b7e47fdbd8cc7acfa 160000
--- a/Tests/Data
+++ b/Tests/Data
@@ -1 +1 @@
-Subproject commit 6de9e6af6d103d65054b2f2b16f6d6b49a53aeaf
+Subproject commit d6ca72279bdc1275b7ccd98b7e47fdbd8cc7acfa
diff --git a/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp b/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..830711966af4e9543285f804ac6eccb2256b1d17
--- /dev/null
+++ b/Tests/ProcessLib/LiquidFlow/TestLiquidFlowMaterialProperties.cpp
@@ -0,0 +1,101 @@
+/**
+ * \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 TestLiquidFlowMaterialProperties.cpp
+ *
+ * Created on August 18, 2016, 4:53 PM
+ */
+
+#include <gtest/gtest.h>
+#include <memory>
+
+#include "TestTools.h"
+
+#include "MeshLib/Mesh.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+
+#include "ProcessLib/Parameter/SpatialPosition.h"
+#include "ProcessLib/LiquidFlow/LiquidFlowMaterialProperties.h"
+
+#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+
+using namespace ProcessLib::LiquidFlow;
+using namespace MaterialLib::Fluid;
+using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType;
+using Matrix = Eigen::MatrixXd;
+
+TEST(ProcessLibLiquidFlow, checkLiquidFlowMaterialProperties)
+{
+ const char xml[] =
+ "<material_property>"
+ " <fluid>"
+ " <density>"
+ " <type>LiquidDensity</type>"
+ " <temperature0> 273.15 </temperature0>"
+ " <p0> 1.e+5 </p0> "
+ " <bulk_modulus> 2.15e+9 </bulk_modulus>"
+ " <beta> 2.0e-4 </beta> "
+ " <rho0>999.8</rho0>"
+ " </density>"
+ " <viscosity>"
+ " <type>Vogels</type>"
+ " <liquid_type>Water </liquid_type>"
+ " </viscosity>"
+ " </fluid>"
+ " <porous_medium>"
+ " <porous_medium id=\"0\">"
+ " <permeability>"
+ " <values>2.e-10 0. 0. 0. 3.e-10 0. 0. 0. 4.0e-10</values>"
+ " </permeability>"
+ " <porosity>"
+ " <type>Constant</type>"
+ " <value> 0.2 </value> "
+ " </porosity>"
+ " <storage>"
+ " <type>Constant</type>"
+ " <value> 1.e-4 </value>"
+ " </storage>"
+ " </porous_medium>"
+ " </porous_medium>"
+ "</material_property>";
+
+ auto const ptree = readXml(xml);
+ BaseLib::ConfigTree conf(ptree, "", BaseLib::ConfigTree::onerror,
+ BaseLib::ConfigTree::onwarning);
+ auto const& sub_config = conf.getConfigSubtree("material_property");
+
+ MeshLib::Properties dummy_property;
+ auto const& dummy_property_vector =
+ dummy_property.createNewPropertyVector<int>(
+ "MaterialIDs", MeshLib::MeshItemType::Cell, 1);
+
+ const bool has_material_ids = false;
+ LiquidFlowMaterialProperties lprop(sub_config, has_material_ids,
+ *dummy_property_vector);
+
+ ProcessLib::SpatialPosition pos;
+ pos.setElementID(0);
+
+ // Check permeability
+ const Eigen::MatrixXd& perm = lprop.getPermeability(0., pos, 1);
+ ASSERT_EQ(2.e-10, perm(0, 0));
+ ASSERT_EQ(0., perm(0, 1));
+ ASSERT_EQ(0., perm(0, 2));
+ ASSERT_EQ(0., perm(1, 0));
+ ASSERT_EQ(3.e-10, perm(1, 1));
+ ASSERT_EQ(0., perm(1, 2));
+ ASSERT_EQ(0., perm(2, 0));
+ ASSERT_EQ(0., perm(2, 1));
+ ASSERT_EQ(4.e-10, perm(2, 2));
+
+ const double T = 273.15 + 60.0;
+ const double p = 1.e+6;
+ const double mass_coef = lprop.getMassCoefficient(0., pos, p, T, 0., 0.);
+ ASSERT_NEAR(0.000100000093, mass_coef, 1.e-10);
+}