diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 72eb86bf8050626414cbb5b382942ad8dae1c9aa..bf53efc64e426a4423472ab7e2d9029773e71f19 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -42,6 +42,7 @@
#include "ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.h"
#include "ProcessLib/TES/CreateTESProcess.h"
#include "ProcessLib/HT/CreateHTProcess.h"
+#include "ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h"
namespace detail
{
@@ -311,8 +312,9 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
else if (type == "LIQUID_FLOW")
{
process = ProcessLib::LiquidFlow::createLiquidFlowProcess(
- *_mesh_vec[0], std::move(jacobian_assembler), _process_variables,
- _parameters, integration_order, process_config);
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config);
}
else if (type == "TES")
{
@@ -408,6 +410,16 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
_process_variables, _parameters, integration_order,
process_config, _curves);
}
+
+ else if (type == "TWOPHASE_FLOW_PP")
+ {
+ process =
+ ProcessLib::TwoPhaseFlowWithPP::CreateTwoPhaseFlowWithPPProcess(
+ *_mesh_vec[0], std::move(jacobian_assembler),
+ _process_variables, _parameters, integration_order,
+ process_config, _curves);
+ }
+
else
{
OGS_FATAL("Unknown process type: %s", type.c_str());
diff --git a/Applications/CLI/Tests.cmake b/Applications/CLI/Tests.cmake
index 2fc03f8b4a3b09b6eb861ef6c1c400e9dfbb437b..18cebaadff4dab757e9dc94577c5f6d3b7fe55ad 100644
--- a/Applications/CLI/Tests.cmake
+++ b/Applications/CLI/Tests.cmake
@@ -363,7 +363,7 @@ if(NOT OGS_USE_MPI)
)
AddTest(
- NAME LARGE_2D_RichardsFlow_h_us_quad
+ NAME LARGE_2D_RichardsFlow_h_us_quad
PATH Parabolic/Richards
EXECUTABLE ogs
EXECUTABLE_ARGS RichardsFlow_2d.prj
@@ -372,6 +372,30 @@ if(NOT OGS_USE_MPI)
DIFF_DATA
h_us_quad_1000.vtu richards_pcs_0_ts_100_t_100.000000.vtu PRESSURE1 pressure
)
+ AddTest(
+ NAME 2D_TwoPhase_PP_Lia_quad
+ PATH Parabolic/TwoPhaseFlowPP/Liakopoulos
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS Twophase_Lia_quad2_short.prj
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-12
+ DIFF_DATA
+ Lia_20.vtu twophaseflow_pcs_0_ts_119_t_20.000000.vtu capillary_pressure capillary_pressure
+ Lia_20.vtu twophaseflow_pcs_0_ts_119_t_20.000000.vtu gas_pressure gas_pressure
+ Lia_20.vtu twophaseflow_pcs_0_ts_119_t_20.000000.vtu saturation saturation
+ )
+ AddTest(
+ NAME LARGE_2D_TwoPhase_PP_Lia_quad
+ PATH Parabolic/TwoPhaseFlowPP/Liakopoulos
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS Twophase_Lia_quad2_large.prj
+ TESTER vtkdiff
+ ABSTOL 1e-8 RELTOL 1e-12
+ DIFF_DATA
+ Lia_1000.vtu twophaseflow_pcs_0_ts_1180_t_1000.000000.vtu capillary_pressure capillary_pressure
+ Lia_1000.vtu twophaseflow_pcs_0_ts_1180_t_1000.000000.vtu gas_pressure gas_pressure
+ Lia_1000.vtu twophaseflow_pcs_0_ts_1180_t_1000.000000.vtu saturation saturation
+ )
AddTest(
NAME LARGE_2D_ThermalConvection_constviscosity
diff --git a/MaterialLib/CMakeLists.txt b/MaterialLib/CMakeLists.txt
index 123a4514098c4718a943572ed5176e64aa4dad7c..9832e12707e8abb532b524eab9f2aef10cc103ab 100644
--- a/MaterialLib/CMakeLists.txt
+++ b/MaterialLib/CMakeLists.txt
@@ -13,6 +13,8 @@ append_source_files(SOURCES PorousMedium/Storage)
append_source_files(SOURCES PorousMedium/Permeability)
append_source_files(SOURCES PorousMedium/UnsaturatedProperty/CapillaryPressure)
+append_source_files(SOURCES TwoPhaseModels)
+
add_library(MaterialLib ${SOURCES} )
target_link_libraries(MaterialLib
BaseLib
diff --git a/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp b/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0e72cac9a809ba5f0dd5f9834a49a5eead0264fb
--- /dev/null
+++ b/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.cpp
@@ -0,0 +1,99 @@
+/**
+ * \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
+ *
+ */
+
+#include "CreateTwoPhaseFlowMaterialProperties.h"
+#include <logog/include/logog.hpp>
+#include "BaseLib/reorderVector.h"
+#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/PropertyVector.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Parameter/SpatialPosition.h"
+#include "TwoPhaseFlowWithPPMaterialProperties.h"
+
+namespace MaterialLib
+{
+namespace TwoPhaseFlowWithPP
+{
+std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>
+CreateTwoPhaseFlowMaterialProperties(
+ BaseLib::ConfigTree const& config,
+ bool const has_material_ids,
+ MeshLib::PropertyVector<int> const& material_ids)
+{
+ DBUG("Reading material properties of two-phase 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__liquid_density}
+ auto const& rho_conf = fluid_config.getConfigSubtree("liquid_density");
+ auto _liquid_density =
+ MaterialLib::Fluid::createFluidDensityModel(rho_conf);
+ //! \ogs_file_param{prj__material_property__gas_density}
+ auto const& rho_gas_conf = fluid_config.getConfigSubtree("gas_density");
+ auto _gas_density =
+ MaterialLib::Fluid::createFluidDensityModel(rho_gas_conf);
+ //! \ogs_file_param{prj__material_property__liquid_viscosity}
+ auto const& mu_conf = fluid_config.getConfigSubtree("liquid_viscosity");
+ auto _viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf);
+ //! \ogs_file_param{prj__material_property__gas_viscosity}
+ auto const& mu_gas_conf = fluid_config.getConfigSubtree("gas_viscosity");
+ auto _gas_viscosity = MaterialLib::Fluid::createViscosityModel(mu_gas_conf);
+
+ // Get porous properties
+ std::vector<int> mat_ids;
+ 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;
+ //! \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& permeability_conf = conf.getConfigSubtree("permeability");
+ _intrinsic_permeability_models.emplace_back(
+ MaterialLib::PorousMedium::createPermeabilityModel(permeability_conf));
+
+ //! \ogs_file_param{prj__material_property__porous_medium__porous_medium__porosity}
+ auto const& porosity_conf = conf.getConfigSubtree("porosity");
+ auto n = MaterialLib::PorousMedium::createPorosityModel(porosity_conf);
+ _porosity_models.emplace_back(std::move(n));
+
+ //! \ogs_file_param{prj__material_property__porous_medium__porous_medium__storage}
+ auto const& storage_conf = conf.getConfigSubtree("storage");
+ auto beta = MaterialLib::PorousMedium::createStorageModel(storage_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);
+
+ return std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>{
+ new TwoPhaseFlowWithPPMaterialProperties{
+ has_material_ids, material_ids, std::move(_liquid_density),
+ std::move(_viscosity), std::move(_gas_density),
+ std::move(_gas_viscosity), _intrinsic_permeability_models,
+ std::move(_porosity_models), std::move(_storage_models)}};
+}
+
+} // end namespace
+} // end namespace
diff --git a/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.h b/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..9e6f31ff7c41d548f8fc1138cf1843e7c24af1e0
--- /dev/null
+++ b/MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.h
@@ -0,0 +1,37 @@
+/**
+ * \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
+ *
+ */
+
+#ifndef OGS_CREATETWOPHASEFLOWMATERIALPROPERTIES_H
+#define OGS_CREATETWOPHASEFLOWMATERIALPROPERTIES_H
+
+#include <memory>
+#include "MaterialLib/Fluid/FluidPropertyHeaders.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h"
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialLib
+{
+namespace TwoPhaseFlowWithPP
+{
+std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>
+CreateTwoPhaseFlowMaterialProperties(
+ BaseLib::ConfigTree const& config,
+ bool const has_material_ids,
+ MeshLib::PropertyVector<int> const& material_ids);
+
+} // end namespace
+} // end namespace
+
+#endif /* CREATETWOPHASEFLOWMATERIALPROPERTIES_H */
diff --git a/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..88ab01f1a021aa2c3f3d8777261efb751a6e8f45
--- /dev/null
+++ b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.cpp
@@ -0,0 +1,133 @@
+/**
+ * \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
+ *
+ */
+
+#include "TwoPhaseFlowWithPPMaterialProperties.h"
+#include <logog/include/logog.hpp>
+#include "BaseLib/reorderVector.h"
+#include "MaterialLib/Fluid/FluidProperty.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "MeshLib/Mesh.h"
+#include "MeshLib/PropertyVector.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Parameter/SpatialPosition.h"
+
+namespace MaterialLib
+{
+namespace TwoPhaseFlowWithPP
+{
+TwoPhaseFlowWithPPMaterialProperties::TwoPhaseFlowWithPPMaterialProperties(
+ bool const has_material_ids,
+ MeshLib::PropertyVector<int> const& material_ids,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ liquid_density,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ viscosity,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ gas_density,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ gas_viscosity,
+ 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)
+ : _has_material_ids(has_material_ids),
+ _liquid_density(std::move(liquid_density)),
+ _viscosity(std::move(viscosity)),
+ _gas_density(std::move(gas_density)),
+ _gas_viscosity(std::move(gas_viscosity)),
+ _material_ids(material_ids),
+ _intrinsic_permeability_models(intrinsic_permeability_models),
+ _porosity_models(std::move(porosity_models)),
+ _storage_models(std::move(storage_models))
+{
+ DBUG("Create material properties for Two-Phase flow with PP model.");
+}
+
+void TwoPhaseFlowWithPPMaterialProperties::setMaterialID(
+ const ProcessLib::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 TwoPhaseFlowWithPPMaterialProperties::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 TwoPhaseFlowWithPPMaterialProperties::getGasDensity(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::pg)] = p;
+ return _gas_density->getValue(vars);
+}
+
+double TwoPhaseFlowWithPPMaterialProperties::getDerivGasDensity(
+ 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::pg)] = p;
+
+ return _gas_density->getdValue(
+ vars, MaterialLib::Fluid::PropertyVariableType::pg);
+}
+double TwoPhaseFlowWithPPMaterialProperties::getLiquidViscosity(
+ 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 TwoPhaseFlowWithPPMaterialProperties::getGasViscosity(
+ 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::pg)] = p;
+ return _gas_viscosity->getValue(vars);
+}
+
+Eigen::MatrixXd const& TwoPhaseFlowWithPPMaterialProperties::getPermeability(
+ const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/,
+ const int /*dim*/) const
+{
+ return _intrinsic_permeability_models[_current_material_id];
+}
+
+double TwoPhaseFlowWithPPMaterialProperties::getPorosity(
+ const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/,
+ const double p, const double T, const double porosity_variable) const
+{
+ const double porosity =
+ _porosity_models[_current_material_id]->getValue(porosity_variable, T);
+
+ return porosity;
+}
+
+} // end of namespace
+} // end of namespace
diff --git a/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..aa3220b9ea4039d90f822011464796da465483ad
--- /dev/null
+++ b/MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h
@@ -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
+ *
+ */
+
+#ifndef OGS_TWOPHASEFLOWWITHPPMATERIALPROPERTIES_H
+#define OGS_TWOPHASEFLOWWITHPPMATERIALPROPERTIES_H
+
+#include <iostream>
+#include <memory>
+#include <vector>
+#include "MaterialLib/Fluid/FluidPropertyHeaders.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+
+// TODO
+// The matereial properties for two phase flow process need to be restructured
+// and moved to a better place.
+namespace ProcessLib
+{
+class SpatialPosition;
+}
+
+namespace MeshLib
+{
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
+namespace MaterialLib
+{
+namespace TwoPhaseFlowWithPP
+{
+class TwoPhaseFlowWithPPMaterialProperties
+{
+public:
+ using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType;
+
+ TwoPhaseFlowWithPPMaterialProperties(
+ bool const has_material_ids,
+ MeshLib::PropertyVector<int> const& material_ids,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ liquid_density,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ viscosity,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ gas_density,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ gas_viscosity,
+ 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);
+
+ void setMaterialID(const ProcessLib::SpatialPosition& pos);
+
+ Eigen::MatrixXd const& getPermeability(
+ const double t,
+ const ProcessLib::SpatialPosition& pos,
+ const int dim) const;
+
+ double getPorosity(const double t, const ProcessLib::SpatialPosition& pos,
+ const double p, const double T,
+ const double porosity_variable) const;
+
+ double getLiquidDensity(const double p, const double T) const;
+ double getGasDensity(const double p, const double T) const;
+ double getGasViscosity(const double p, const double T) const;
+ double getLiquidViscosity(const double p, const double T) const;
+ double getDerivGasDensity(double const p, double const T) const;
+
+protected:
+ /// A flag to indicate whether the reference member, _material_ids,
+ /// is not assigned.
+ const bool _has_material_ids;
+
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_density;
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty> _viscosity;
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_density;
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty> _gas_viscosity;
+
+ /** Use two phase 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;
+};
+
+} // end of namespace
+} // end of namespace
+#endif /* TWOPHASEFLOWWITHPPMATERIALPROPERTIES_H */
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index 592ba502580ea8d608bdb92de4a74602bc7cdccd..6425192f1cdf6b7982de28937f4a5784ab91126f 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -17,6 +17,8 @@ APPEND_SOURCE_FILES(SOURCES GroundwaterFlow)
APPEND_SOURCE_FILES(SOURCES LiquidFlow)
APPEND_SOURCE_FILES(SOURCES HydroMechanics)
+APPEND_SOURCE_FILES(SOURCES TwoPhaseFlowWithPP)
+
APPEND_SOURCE_FILES(SOURCES SmallDeformation)
APPEND_SOURCE_FILES(SOURCES LIE/BoundaryCondition)
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..650278df43db14480c28602eb1317fd86641153c
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp
@@ -0,0 +1,134 @@
+/**
+ * \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
+ *
+ */
+#include "CreateTwoPhaseFlowWithPPProcess.h"
+#include <cassert>
+
+#include "MaterialLib/TwoPhaseModels/CreateTwoPhaseFlowMaterialProperties.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+#include "ProcessLib/Parameter/ConstantParameter.h"
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+
+#include "TwoPhaseFlowWithPPProcess.h"
+#include "TwoPhaseFlowWithPPProcessData.h"
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPP
+{
+std::unique_ptr<Process> CreateTwoPhaseFlowWithPPProcess(
+ 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,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves)
+{
+ //! \ogs_file_param{process__type}
+ config.checkConfigParameter("type", "TWOPHASE_FLOW_PP");
+
+ DBUG("Create TwoPhaseFlowProcess with PP model.");
+
+ // Process variable.
+ auto process_variables = findProcessVariables(
+ variables, config,
+ {//! \ogs_file_param_special{process__TWOPHASE_FLOW_PP__process_variables__process_variable}
+ "gas_pressure", "capillary_pressure"});
+
+ SecondaryVariableCollection secondary_variables;
+
+ NumLib::NamedFunctionCaller named_function_caller(
+ {"TwoPhaseFlow_pressure"});
+
+ ProcessLib::parseSecondaryVariables(config, secondary_variables,
+ named_function_caller);
+ // Specific body force
+ Eigen::VectorXd specific_body_force;
+
+ std::vector<double> const b =
+ //! \ogs_file_param_special{process__TWOPHASE_FLOW_PP__specific_body_force}
+ config.getConfigParameter<std::vector<double>>("specific_body_force");
+ assert(b.size() > 0 && b.size() < 4);
+ specific_body_force.resize(b.size());
+ bool const has_gravity = MathLib::toVector(b).norm() > 0;
+ if (has_gravity)
+ std::copy_n(b.data(), b.size(), specific_body_force.data());
+
+ // //! \ogs_file_param{process__TWOPHASE_FLOW_PP__mass_lumping}
+ auto mass_lump = config.getConfigParameter<bool>("mass_lumping");
+
+ //! \ogs_file_param{process__TWOPHASE_FLOW_PP__material_property}
+ auto const& mat_config = config.getConfigSubtree("material_property");
+
+ auto const& mat_ids =
+ mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+
+ std::unique_ptr<
+ MaterialLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPMaterialProperties>
+ material = nullptr;
+
+ if (mat_ids)
+ {
+ INFO("The twophase flow is in heterogeneous porous media.");
+ const bool has_material_ids = true;
+ material = MaterialLib::TwoPhaseFlowWithPP::
+ CreateTwoPhaseFlowMaterialProperties(mat_config, has_material_ids,
+ *mat_ids);
+ TwoPhaseFlowWithPPProcessData process_data{
+ specific_body_force,
+ has_gravity,
+ mass_lump,
+ std::move(material),
+ *curves.at("curve_PC_S"),
+ *curves.at("curve_S_Krel_wet"),
+ *curves.at("curve_S_Krel_nonwet")};
+ return std::unique_ptr<Process>{new TwoPhaseFlowWithPPProcess{
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(process_data),
+ std::move(secondary_variables), std::move(named_function_caller),
+ mat_config, curves}};
+ }
+ else
+ {
+ INFO("The twophase flow is in homogeneous porous media.");
+
+ MeshLib::Properties dummy_property;
+
+ 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;
+ material = MaterialLib::TwoPhaseFlowWithPP::
+ CreateTwoPhaseFlowMaterialProperties(mat_config, has_material_ids,
+ *dummy_property_vector);
+ TwoPhaseFlowWithPPProcessData process_data{
+ specific_body_force,
+ has_gravity,
+ mass_lump,
+ std::move(material),
+ *curves.at("curve_PC_S"),
+ *curves.at("curve_S_Krel_wet"),
+ *curves.at("curve_S_Krel_nonwet")};
+ return std::unique_ptr<Process>{new TwoPhaseFlowWithPPProcess{
+ mesh, std::move(jacobian_assembler), parameters, integration_order,
+ std::move(process_variables), std::move(process_data),
+ std::move(secondary_variables), std::move(named_function_caller),
+ mat_config, curves}};
+ }
+}
+
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..d7e14e7c496f728085ce4396bb47523863518819
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.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
+ *
+ */
+
+#ifndef OGS_CREATETWOPHASEFLOWWITHPPPROCESS_H
+#define OGS_CREATETWOPHASEFLOWWITHPPPROCESS_H
+
+#include <memory>
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPP
+{
+std::unique_ptr<Process> CreateTwoPhaseFlowWithPPProcess(
+ 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,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves);
+} // end of namespace
+} // end of namespace
+
+#endif /* CREATETWOPHASEFLOWWITHPPPROCESS_H */
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..a47d27f853bf9e566f0fa7d8c74086d16bf4e8cb
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
@@ -0,0 +1,208 @@
+/**
+ * \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
+ *
+ */
+
+#ifndef OGS_TWOPHASEFLOWWITHPPLOCALASSEMBLER_IMPL_H
+#define OGS_TWOPHASEFLOWWITHPPLOCALASSEMBLER_IMPL_H
+
+#include "TwoPhaseFlowWithPPLocalAssembler.h"
+
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "NumLib/Function/Interpolation.h"
+#include "TwoPhaseFlowWithPPProcessData.h"
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPP
+{
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void TwoPhaseFlowWithPPLocalAssembler<
+ 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)
+{
+ auto const local_matrix_size = local_x.size();
+
+ assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
+
+ auto local_M = MathLib::createZeroedMatrix<LocalMatrixType>(
+ local_M_data, local_matrix_size, local_matrix_size);
+ auto local_K = MathLib::createZeroedMatrix<LocalMatrixType>(
+ local_K_data, local_matrix_size, local_matrix_size);
+ auto local_b = MathLib::createZeroedVector<LocalVectorType>(
+ local_b_data, local_matrix_size);
+
+ NodalMatrixType mass_operator;
+ mass_operator.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
+
+ auto Mgp =
+ local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
+ nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
+ auto Mgpc = local_M.template block<nonwet_pressure_size, cap_pressure_size>(
+ nonwet_pressure_matrix_index, cap_pressure_matrix_index);
+
+ auto Mlp = local_M.template block<cap_pressure_size, nonwet_pressure_size>(
+ cap_pressure_matrix_index, nonwet_pressure_matrix_index);
+
+ auto Mlpc = local_M.template block<cap_pressure_size, cap_pressure_size>(
+ cap_pressure_matrix_index, cap_pressure_matrix_index);
+
+ NodalMatrixType laplace_operator;
+ laplace_operator.setZero(ShapeFunction::NPOINTS, ShapeFunction::NPOINTS);
+
+ auto Kgp =
+ local_K.template block<nonwet_pressure_size, nonwet_pressure_size>(
+ nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
+
+ auto Kgpc = local_K.template block<nonwet_pressure_size, cap_pressure_size>(
+ nonwet_pressure_matrix_index, cap_pressure_matrix_index);
+
+ auto Klp = local_K.template block<cap_pressure_size, nonwet_pressure_size>(
+ cap_pressure_matrix_index, nonwet_pressure_matrix_index);
+
+ auto Klpc = local_K.template block<cap_pressure_size, cap_pressure_size>(
+ cap_pressure_matrix_index, cap_pressure_matrix_index);
+
+ auto Bg = local_b.template segment<nonwet_pressure_size>(
+ nonwet_pressure_matrix_index);
+
+ auto Bl =
+ local_b.template segment<cap_pressure_size>(cap_pressure_matrix_index);
+
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+
+ SpatialPosition pos;
+ pos.setElementID(_element.getID());
+ _process_data._material->setMaterialID(pos);
+
+ const Eigen::MatrixXd& perm = _process_data._material->getPermeability(
+ t, pos, _element.getDimension());
+ assert(perm.rows() == GlobalDim || perm.rows() == 1);
+ GlobalDimMatrixType permeability =
+ GlobalDimMatrixType::Zero(GlobalDim, GlobalDim);
+ if (perm.rows() == GlobalDim)
+ permeability = perm;
+ else if (perm.rows() == 1)
+ permeability.diagonal().setConstant(perm(0, 0));
+ MathLib::PiecewiseLinearInterpolation const& interpolated_Pc =
+ _process_data._interpolated_Pc;
+ MathLib::PiecewiseLinearInterpolation const& interpolated_Kr_wet =
+ _process_data._interpolated_Kr_wet;
+ MathLib::PiecewiseLinearInterpolation const& interpolated_Kr_nonwet =
+ _process_data._interpolated_Kr_nonwet;
+
+ // Note: currently only isothermal case is considered, so the temperature is
+ // assumed to be const
+ // the variation of temperature will be taken into account in future
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& sm = _shape_matrices[ip];
+
+ double pc_int_pt = 0.;
+ double pg_int_pt = 0.;
+ NumLib::shapeFunctionInterpolate(local_x, sm.N, pg_int_pt, pc_int_pt);
+
+ _pressure_wetting[ip] = pg_int_pt - pc_int_pt;
+
+ auto const& wp = _integration_method.getWeightedPoint(ip);
+ const double integration_factor =
+ sm.integralMeasure * sm.detJ * wp.getWeight();
+
+ double const rho_gas =
+ _process_data._material->getGasDensity(pg_int_pt, _temperature);
+ double const rho_w = _process_data._material->getLiquidDensity(
+ _pressure_wetting[ip], _temperature);
+
+ double const Sw =
+ (pc_int_pt < 0) ? 1.0 : interpolated_Pc.getValue(pc_int_pt);
+
+ _saturation[ip] = Sw;
+ double dSwdPc = interpolated_Pc.getDerivative(pc_int_pt);
+ if (pc_int_pt > interpolated_Pc.getSupportMax())
+ dSwdPc =
+ interpolated_Pc.getDerivative(interpolated_Pc.getSupportMax());
+ else if (pc_int_pt < interpolated_Pc.getSupportMin())
+ dSwdPc =
+ interpolated_Pc.getDerivative(interpolated_Pc.getSupportMin());
+
+ double const porosity = _process_data._material->getPorosity(
+ t, pos, pg_int_pt, _temperature, 0);
+
+ // Assemble M matrix
+ // nonwetting
+ double const drhogas_dpg = _process_data._material->getDerivGasDensity(
+ pg_int_pt, _temperature);
+
+ mass_operator.noalias() = sm.N.transpose() * sm.N * integration_factor;
+
+ Mgp.noalias() += porosity * (1 - Sw) * drhogas_dpg * mass_operator;
+ Mgpc.noalias() += -porosity * rho_gas * dSwdPc * mass_operator;
+ Mlp.noalias() += 0.0 * mass_operator;
+ Mlpc.noalias() += porosity * dSwdPc * rho_w * mass_operator;
+
+ // Assemble M matrix
+ // nonwet
+ double const k_rel_G = interpolated_Kr_nonwet.getValue(Sw);
+ double const mu_gas =
+ _process_data._material->getGasViscosity(pg_int_pt, _temperature);
+ double const lambda_G = k_rel_G / mu_gas;
+
+ // wet
+ double const k_rel_L = interpolated_Kr_wet.getValue(Sw);
+ double const mu_liquid = _process_data._material->getLiquidViscosity(
+ _pressure_wetting[ip], _temperature);
+ double const lambda_L = k_rel_L / mu_liquid;
+
+ laplace_operator.noalias() =
+ sm.dNdx.transpose() * permeability * sm.dNdx * integration_factor;
+
+ Kgp.noalias() += rho_gas * lambda_G * laplace_operator;
+ Kgpc.noalias() += 0.0 * laplace_operator;
+ Klp.noalias() += rho_w * lambda_L * laplace_operator;
+ Klpc.noalias() += -rho_w * lambda_L * laplace_operator;
+
+ if (_process_data._has_gravity)
+ {
+ auto const& b = _process_data._specific_body_force;
+ Bg.noalias() += rho_gas * rho_gas * lambda_G * sm.dNdx.transpose() *
+ permeability * b * integration_factor;
+ Bl.noalias() += rho_w * rho_w * lambda_L * sm.dNdx.transpose() *
+ permeability * b * integration_factor;
+
+ } // end of has gravity
+ } // end of GP
+ if (_process_data._has_mass_lumping)
+ {
+ for (unsigned row = 0; row < Mgpc.cols(); row++)
+ {
+ for (unsigned column = 0; column < Mgpc.cols(); column++)
+ {
+ if (row != column)
+ {
+ Mgpc(row, row) += Mgpc(row, column);
+ Mgpc(row, column) = 0.0;
+ Mgp(row, row) += Mgp(row, column);
+ Mgp(row, column) = 0.0;
+ Mlpc(row, row) += Mlpc(row, column);
+ Mlpc(row, column) = 0.0;
+ Mlp(row, row) += Mlp(row, column);
+ Mlp(row, column) = 0.0;
+ }
+ }
+ }
+ } // end of mass-lumping
+}
+
+} // end of namespace
+} // end of namespace
+
+#endif
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..136fa3a6d6c3770234bdbb40120cba021976861d
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
@@ -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
+ *
+ */
+
+#ifndef OGS_TWOPHASEFLOWWITHPPLOCALASSEMBLER_H
+#define OGS_TWOPHASEFLOWWITHPPLOCALASSEMBLER_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 "MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h"
+#include "TwoPhaseFlowWithPPProcessData.h"
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPP
+{
+const unsigned NUM_NODAL_DOF = 2;
+
+class TwoPhaseFlowWithPPLocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+public:
+ virtual std::vector<double> const& getIntPtSaturation(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtWettingPressure(
+ std::vector<double>& /*cache*/) const = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class TwoPhaseFlowWithPPLocalAssembler
+ : public TwoPhaseFlowWithPPLocalAssemblerInterface
+{
+ using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;
+ using ShapeMatrices = typename ShapeMatricesType::ShapeMatrices;
+
+ using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits<
+ ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
+
+ using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
+ using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
+ using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
+ using GlobalDimVectorType = typename ShapeMatricesType::GlobalDimVectorType;
+ using LocalMatrixType = typename LocalAssemblerTraits::LocalMatrix;
+ using LocalVectorType = typename LocalAssemblerTraits::LocalVector;
+
+public:
+ TwoPhaseFlowWithPPLocalAssembler(
+ MeshLib::Element const& element,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ TwoPhaseFlowWithPPProcessData const& process_data)
+ : _element(element),
+ _integration_method(integration_order),
+ _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ element, is_axially_symmetric, _integration_method)),
+ _process_data(process_data),
+ _saturation(
+ std::vector<double>(_integration_method.getNumberOfPoints())),
+ _pressure_wetting(
+ std::vector<double>(_integration_method.getNumberOfPoints()))
+ {
+ }
+
+ 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& getIntPtSaturation(
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_saturation.size() > 0);
+ return _saturation;
+ }
+
+ std::vector<double> const& getIntPtWettingPressure(
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_pressure_wetting.size() > 0);
+ return _pressure_wetting;
+ }
+
+private:
+ MeshLib::Element const& _element;
+
+ IntegrationMethod const _integration_method;
+ std::vector<ShapeMatrices> _shape_matrices;
+
+ TwoPhaseFlowWithPPProcessData const& _process_data;
+
+ // Note: currently only isothermal case is considered, so the temperature is
+ // assumed to be const
+ // the variation of temperature will be taken into account in future
+ double _temperature = 293.15;
+ std::vector<double> _saturation;
+ std::vector<double> _pressure_wetting;
+ static const int nonwet_pressure_coeff_index = 0;
+ static const int cap_pressure_coeff_index = 1;
+
+ static const int nonwet_pressure_matrix_index = 0;
+ static const int cap_pressure_matrix_index = ShapeFunction::NPOINTS;
+
+ static const int nonwet_pressure_size = ShapeFunction::NPOINTS;
+ static const int cap_pressure_size = ShapeFunction::NPOINTS;
+};
+
+} // end of namespace
+} // end of namespace
+
+#include "TwoPhaseFlowWithPPLocalAssembler-impl.h"
+
+#endif /* TWOPHASEFLOWWITHPPLOCALASSEMBLER_H */
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0f1a7cc7a3949d474bb75bff5cb68416d11e6e25
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -0,0 +1,99 @@
+/**
+ * \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
+ *
+ */
+
+#include "TwoPhaseFlowWithPPProcess.h"
+
+#include <cassert>
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "MeshLib/PropertyVector.h"
+
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "TwoPhaseFlowWithPPLocalAssembler.h"
+
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPP
+{
+TwoPhaseFlowWithPPProcess::TwoPhaseFlowWithPPProcess(
+ 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,
+ TwoPhaseFlowWithPPProcessData&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller,
+ BaseLib::ConfigTree const& config,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves)
+ : Process(mesh, std::move(jacobian_assembler), parameters,
+ integration_order, std::move(process_variables),
+ std::move(secondary_variables), std::move(named_function_caller)),
+ _process_data(std::move(process_data))
+{
+ DBUG("Create TwoPhaseFlowProcess with PP model.");
+}
+
+void TwoPhaseFlowWithPPProcess::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ ProcessLib::ProcessVariable const& pv = getProcessVariables()[0];
+ ProcessLib::createLocalAssemblers<TwoPhaseFlowWithPPLocalAssembler>(
+ mesh.getDimension(), mesh.getElements(), dof_table,
+ pv.getShapeFunctionOrder(), _local_assemblers,
+ mesh.isAxiallySymmetric(), integration_order, _process_data);
+
+ _secondary_variables.addSecondaryVariable(
+ "saturation", 1,
+ makeExtrapolator(
+ getExtrapolator(), _local_assemblers,
+ &TwoPhaseFlowWithPPLocalAssemblerInterface::getIntPtSaturation));
+
+ _secondary_variables.addSecondaryVariable(
+ "pressure_wetting", 1,
+ makeExtrapolator(getExtrapolator(), _local_assemblers,
+ &TwoPhaseFlowWithPPLocalAssemblerInterface::
+ getIntPtWettingPressure));
+}
+
+void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t,
+ GlobalVector const& x,
+ GlobalMatrix& M,
+ GlobalMatrix& K,
+ GlobalVector& b)
+{
+ DBUG("Assemble TwoPhaseFlowWithPPProcess.");
+ // 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 TwoPhaseFlowWithPPProcess::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 TwoPhaseFlowWithPPProcess.");
+
+ // 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/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..f91c34bc21c5390ba0fb6747e90153a23b73e7f5
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
@@ -0,0 +1,77 @@
+/**
+ * \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
+ *
+ */
+
+#ifndef OGS_TWOPHASEFLOWWITHPPPROCESS_H
+#define OGS_TWOPHASEFLOWWITHPPPROCESS_H
+
+#include "MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h"
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/Process.h"
+#include "TwoPhaseFlowWithPPLocalAssembler.h"
+
+namespace MeshLib
+{
+class Element;
+class Mesh;
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPP
+{
+/**
+ * \brief A class to simulate the two-phase flow process with P-P model in
+ * porous media
+ */
+class TwoPhaseFlowWithPPProcess final : public Process
+{
+public:
+ TwoPhaseFlowWithPPProcess(
+ 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,
+ TwoPhaseFlowWithPPProcessData&& process_data,
+ SecondaryVariableCollection&& secondary_variables,
+ NumLib::NamedFunctionCaller&& named_function_caller,
+ BaseLib::ConfigTree const& config,
+ std::map<std::string,
+ std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
+ curves);
+
+ bool isLinear() const override { return false; }
+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;
+
+ TwoPhaseFlowWithPPProcessData _process_data;
+
+ std::vector<std::unique_ptr<TwoPhaseFlowWithPPLocalAssemblerInterface>>
+ _local_assemblers;
+};
+
+} // end of namespace
+} // end of namespace
+
+#endif /* TWOPHASEFLOWWITHPPPROCESS_H */
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h
new file mode 100644
index 0000000000000000000000000000000000000000..ea612754122ba90487f5982275dfe7a2d1cfa3dc
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h
@@ -0,0 +1,81 @@
+/**
+ * \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
+ *
+ */
+
+#ifndef PROCESSLIB_TWOPHASEFLOWWITHPP_TWOPHASEFLOWWITHPPPROCESSDATA_H
+#define PROCESSLIB_TWOPHASEFLOWWITHPP_TWOPHASEFLOWWITHPPPROCESSDATA_H
+
+namespace MeshLib
+{
+class Element;
+}
+
+namespace ProcessLib
+{
+template <typename T>
+struct Parameter;
+
+namespace TwoPhaseFlowWithPP
+{
+struct TwoPhaseFlowWithPPProcessData
+{
+ TwoPhaseFlowWithPPProcessData(
+ Eigen::VectorXd const specific_body_force_,
+ bool const has_gravity_,
+ bool const has_mass_lumping_,
+ std::unique_ptr<MaterialLib::TwoPhaseFlowWithPP::
+ TwoPhaseFlowWithPPMaterialProperties>&& material_,
+ MathLib::PiecewiseLinearInterpolation const& interpolated_Pc_,
+ MathLib::PiecewiseLinearInterpolation const& interpolated_Kr_wet_,
+ MathLib::PiecewiseLinearInterpolation const& interpolated_Kr_nonwet_)
+ : _specific_body_force(specific_body_force_),
+ _has_gravity(has_gravity_),
+ _has_mass_lumping(has_mass_lumping_),
+ _material(std::move(material_)),
+ _interpolated_Pc(interpolated_Pc_),
+ _interpolated_Kr_wet(interpolated_Kr_wet_),
+ _interpolated_Kr_nonwet(interpolated_Kr_nonwet_)
+ {
+ }
+
+ TwoPhaseFlowWithPPProcessData(TwoPhaseFlowWithPPProcessData&& other)
+ : _specific_body_force(other._specific_body_force),
+ _has_gravity(other._has_gravity),
+ _has_mass_lumping(other._has_mass_lumping),
+ _material(std::move(other._material)),
+ _interpolated_Pc(other._interpolated_Pc),
+ _interpolated_Kr_wet(other._interpolated_Kr_wet),
+ _interpolated_Kr_nonwet(other._interpolated_Kr_nonwet)
+ {
+ }
+
+ //! Copies are forbidden.
+ TwoPhaseFlowWithPPProcessData(TwoPhaseFlowWithPPProcessData const&) =
+ delete;
+
+ //! Assignments are not needed.
+ void operator=(TwoPhaseFlowWithPPProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(TwoPhaseFlowWithPPProcessData&&) = delete;
+ Eigen::VectorXd const _specific_body_force;
+
+ bool const _has_gravity;
+ bool const _has_mass_lumping;
+ std::unique_ptr<
+ MaterialLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPMaterialProperties>
+ _material;
+ MathLib::PiecewiseLinearInterpolation const& _interpolated_Pc;
+ MathLib::PiecewiseLinearInterpolation const& _interpolated_Kr_wet;
+ MathLib::PiecewiseLinearInterpolation const& _interpolated_Kr_nonwet;
+};
+
+} // namespace TwoPhaseFlowWithPP
+} // namespace ProcessLib
+
+#endif // PROCESSLIB_TwoPhaseFlowWithPP_TwoPhaseFlowWithPPPROCESSDATA_H