diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 2436b1a083b93c709cd2d1e725ae511e0c878735..d3d46f9b44269b2af7d50493a89857026d2a860f 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -437,7 +437,7 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
else if (type == "TWOPHASE_FLOW_PP")
{
process =
- ProcessLib::TwoPhaseFlowWithPP::CreateTwoPhaseFlowWithPPProcess(
+ ProcessLib::TwoPhaseFlowWithPP::createTwoPhaseFlowWithPPProcess(
*_mesh_vec[0], std::move(jacobian_assembler),
_process_variables, _parameters, integration_order,
process_config, _curves);
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/i_relative_permeability.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/i_relative_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..e3e0c731598545449b3007e7346da0d21e57e6fd
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/i_relative_permeability.md
@@ -0,0 +1 @@
+A tag for the relative permeability model, which holds two ids: one is for the Nonwetting phase and the other is for the Wetting phase.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/a_id.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/a_id.md
new file mode 100644
index 0000000000000000000000000000000000000000..a64a7806df41e0c3e78448f823887d29b49859c8
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/a_id.md
@@ -0,0 +1 @@
+A tag for the relative permeability model id: 0 is Nonwetting phase and 1 is Wetting phase.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/i_relative_permeability.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/i_relative_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..ea0efd1bb5564eb77343a383d8314799cd0443b1
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/i_relative_permeability.md
@@ -0,0 +1 @@
+Defines the relative permeability model.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/t_capillary_pressure.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/t_capillary_pressure.md
new file mode 100644
index 0000000000000000000000000000000000000000..19d854870d0c3537197cc8edb349b03503fcbede
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/porous_medium/porous_medium/t_capillary_pressure.md
@@ -0,0 +1 @@
+Capillary pressure model.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_density.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_density.md
new file mode 100644
index 0000000000000000000000000000000000000000..bcf97b55a9e2bf33941d089ec03b97f90d562b1f
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_density.md
@@ -0,0 +1 @@
+Density model of the gas phase.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_viscosity.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_viscosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..a02e98700f601132ce1e4e545b1f87d2d0eb5388
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_gas_viscosity.md
@@ -0,0 +1 @@
+Viscosity model of the gas phase.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_density.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_density.md
new file mode 100644
index 0000000000000000000000000000000000000000..5c8db5dec3d6f9cc8f530470cbca94b43e23411c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_density.md
@@ -0,0 +1 @@
+Density model of the liquid phase.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_viscosity.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_viscosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..0451570185b81e26ecab7dad8eb9b7af93776891
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/material_property/t_liquid_viscosity.md
@@ -0,0 +1 @@
+Viscosity model of the liquid phase.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_mass_lumping.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_mass_lumping.md
index 0d6032b0aad3bc8e76544c5c6669aacc6d24ef6b..2cfa81d1102216d98163c46246a1c925df4d43f2 100644
--- a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_mass_lumping.md
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_mass_lumping.md
@@ -1 +1 @@
-\copydoc ProcessLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPProcessData::_has_mass_lumping
+\copydoc ProcessLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPProcessData::has_mass_lumping
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_specific_body_force.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_specific_body_force.md
index 7437fbf0be7c14f306021d29126880811985c36c..c1658dbc4ddc31d43e17c96ebc717ed667bac0ea 100644
--- a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_specific_body_force.md
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_specific_body_force.md
@@ -1 +1 @@
-\copydoc ProcessLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPProcessData::_specific_body_force
+\copydoc ProcessLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPProcessData::specific_body_force
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_temperature.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_temperature.md
new file mode 100644
index 0000000000000000000000000000000000000000..ee031866715fe568cec9def41224bef991df9201
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PP/t_temperature.md
@@ -0,0 +1,2 @@
+An input of a reference temperature (K) for two-phase flow process.
+It is fixed to be constant for the simulation due to an isothermal assumption.
diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
index 4fec226af6f7e830d37b4d312d4a74ac95b4436e..df5ac0e9b9a69419ce52f501ed5366496215e580 100644
--- a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
+++ b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
@@ -95,7 +95,7 @@ double PiecewiseLinearInterpolation::getDerivative(
interval_idx = 1;
}
- if (interval_idx > 2 && interval_idx < _supp_pnts.size() - 1)
+ if (interval_idx > 1 && interval_idx < _supp_pnts.size() - 2)
{
// left and right support points.
double const x_ll = _supp_pnts[interval_idx - 2];
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CMakeLists.txt b/ProcessLib/TwoPhaseFlowWithPP/CMakeLists.txt
index b21e3b6008d9bd837cc2964ad2d99a1da90f3fa1..12a9b6994522f0f4cdbac7d00cdb3d08653d3369 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/CMakeLists.txt
+++ b/ProcessLib/TwoPhaseFlowWithPP/CMakeLists.txt
@@ -1,26 +1,69 @@
AddTest(
- NAME 2D_TwoPhase_PP_Lia_quad
+ NAME 2D_TwoPhase_PP_Lia_quad_1
PATH Parabolic/TwoPhaseFlowPP/Liakopoulos
EXECUTABLE ogs
- EXECUTABLE_ARGS Twophase_Lia_quad2_short.prj
+ EXECUTABLE_ARGS TwoPhase_Lia_quad_short.prj
TESTER vtkdiff
REQUIREMENTS NOT OGS_USE_MPI
- ABSTOL 1e-8 RELTOL 1e-12
+ ABSTOL 1e-2 RELTOL 1e-4
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
+ h2_Liako_20.vtu twophaseflow_pcs_0_ts_218_t_20.000000.vtu saturation saturation
)
AddTest(
- NAME LARGE_2D_TwoPhase_PP_Lia_quad
+ NAME 2D_TwoPhase_PP_Lia_quad_2
PATH Parabolic/TwoPhaseFlowPP/Liakopoulos
EXECUTABLE ogs
- EXECUTABLE_ARGS Twophase_Lia_quad2_large.prj
+ EXECUTABLE_ARGS TwoPhase_Lia_quad_short.prj
TESTER vtkdiff
REQUIREMENTS NOT OGS_USE_MPI
- ABSTOL 1e-8 RELTOL 1e-12
+ ABSTOL 20 RELTOL 1e-3
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
+ h2_Liako_20.vtu twophaseflow_pcs_0_ts_218_t_20.000000.vtu capillary_pressure capillary_pressure
+ h2_Liako_20.vtu twophaseflow_pcs_0_ts_218_t_20.000000.vtu gas_pressure gas_pressure
+)
+AddTest(
+ NAME LARGE_2D_TwoPhase_PP_Lia_quad_1
+ PATH Parabolic/TwoPhaseFlowPP/Liakopoulos
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS TwoPhase_Lia_quad_large.prj
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ ABSTOL 1e-2 RELTOL 1e-3
+ DIFF_DATA
+ h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu SATURATION1 saturation
+)
+AddTest(
+ NAME LARGE_2D_TwoPhase_PP_Lia_quad_2
+ PATH Parabolic/TwoPhaseFlowPP/Liakopoulos
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS TwoPhase_Lia_quad_large.prj
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ ABSTOL 20 RELTOL 1e-2
+ DIFF_DATA
+ h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu PRESSURE1 capillary_pressure
+ h2_Liako_1198.vtu twophaseflow_pcs_0_ts_1198_t_1000.000000.vtu PRESSURE2 gas_pressure
+)
+AddTest(
+ NAME 1D_TwoPhase_PP_mcwt_1
+ PATH Parabolic/TwoPhaseFlowPP/McWorts
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS TwoPhase_mcwt_line.prj
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ ABSTOL 1e-2 RELTOL 1e-4
+ DIFF_DATA
+ mcwt_1000.vtu twophaseflow_pcs_0_ts_519_t_1000.000000.vtu SATURATION1 saturation
+)
+AddTest(
+ NAME 1D_TwoPhase_PP_mcwt_2
+ PATH Parabolic/TwoPhaseFlowPP/McWorts
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS TwoPhase_mcwt_line.prj
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ ABSTOL 10 RELTOL 1e-3
+ DIFF_DATA
+ mcwt_1000.vtu twophaseflow_pcs_0_ts_519_t_1000.000000.vtu PRESSURE1 capillary_pressure
+ mcwt_1000.vtu twophaseflow_pcs_0_ts_519_t_1000.000000.vtu PRESSURE2 gas_pressure
)
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5c5f353d016a0a5c0f96acc098a320c6db0467b1
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.cpp
@@ -0,0 +1,139 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, 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 "CreateTwoPhaseFlowWithPPMaterialProperties.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 "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.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 ProcessLib
+{
+namespace TwoPhaseFlowWithPP
+{
+std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>
+createTwoPhaseFlowWithPPMaterialProperties(
+ BaseLib::ConfigTree const& config,
+ boost::optional<MeshLib::PropertyVector<int> const&>
+ material_ids)
+{
+ DBUG("Reading material properties of two-phase flow process.");
+
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__fluid}
+ auto const& fluid_config = config.getConfigSubtree("fluid");
+
+ // Get fluid properties
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__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__processes__process__TWOPHASE_FLOW_PP__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__processes__process__TWOPHASE_FLOW_PP__material_property__liquid_viscosity}
+ auto const& mu_conf = fluid_config.getConfigSubtree("liquid_viscosity");
+ auto liquid_viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf);
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__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<int> mat_krel_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;
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>>
+ capillary_pressure_models;
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>
+ relative_permeability_models;
+
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium}
+ auto const& poro_config = config.getConfigSubtree("porous_medium");
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium}
+ for (auto const& conf : poro_config.getConfigSubtreeList("porous_medium"))
+ {
+ //! \ogs_file_attr{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__id}
+ auto const id = conf.getConfigAttributeOptional<int>("id");
+ mat_ids.push_back(*id);
+
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__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__processes__process__TWOPHASE_FLOW_PP__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__processes__process__TWOPHASE_FLOW_PP__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));
+
+ auto const& capillary_pressure_conf =
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__capillary_pressure}
+ conf.getConfigSubtree("capillary_pressure");
+ auto pc = MaterialLib::PorousMedium::createCapillaryPressureModel(
+ capillary_pressure_conf);
+ capillary_pressure_models.emplace_back(std::move(pc));
+
+ auto const& krel_config =
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__relative_permeability}
+ conf.getConfigSubtree("relative_permeability");
+ for (
+ auto const& krel_conf :
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__relative_permeability__relative_permeability}
+ krel_config.getConfigSubtreeList("relative_permeability"))
+ {
+ auto const krel_id =
+ //! \ogs_file_attr{prj__processes__process__TWOPHASE_FLOW_PP__material_property__porous_medium__porous_medium__relative_permeability__relative_permeability__id}
+ krel_conf.getConfigAttributeOptional<int>("id");
+ mat_krel_ids.push_back(*krel_id);
+ auto krel_n =
+ MaterialLib::PorousMedium::createRelativePermeabilityModel(
+ krel_conf);
+ relative_permeability_models.emplace_back(std::move(krel_n));
+ }
+ BaseLib::reorderVector(relative_permeability_models, mat_krel_ids);
+ }
+
+ 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{
+ material_ids, std::move(liquid_density),
+ std::move(liquid_viscosity), std::move(gas_density),
+ std::move(gas_viscosity), intrinsic_permeability_models,
+ std::move(porosity_models), std::move(storage_models),
+ std::move(capillary_pressure_models),
+ std::move(relative_permeability_models)}};
+}
+
+} // end namespace
+} // end namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.h b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..a00d066f3821d76ed6be1b24b27bb04c612a01e3
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPMaterialProperties.h
@@ -0,0 +1,30 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+
+#include "TwoPhaseFlowWithPPMaterialProperties.h"
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPP
+{
+std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>
+createTwoPhaseFlowWithPPMaterialProperties(
+ BaseLib::ConfigTree const& config,
+ boost::optional<MeshLib::PropertyVector<int> const&>
+ material_ids);
+
+} // end namespace
+} // end namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp
index 460edafca73b46490d4048747991d0ed721be349..03375787171a34c465f5c70acbbc1b8131f00a0e 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.cpp
@@ -9,20 +9,20 @@
#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 "CreateTwoPhaseFlowWithPPMaterialProperties.h"
+#include "TwoPhaseFlowWithPPMaterialProperties.h"
#include "TwoPhaseFlowWithPPProcess.h"
#include "TwoPhaseFlowWithPPProcessData.h"
-
namespace ProcessLib
{
namespace TwoPhaseFlowWithPP
{
-std::unique_ptr<Process> CreateTwoPhaseFlowWithPPProcess(
+std::unique_ptr<Process> createTwoPhaseFlowWithPPProcess(
MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
@@ -37,9 +37,6 @@ std::unique_ptr<Process> CreateTwoPhaseFlowWithPPProcess(
config.checkConfigParameter("type", "TWOPHASE_FLOW_PP");
DBUG("Create TwoPhaseFlowProcess with PP model.");
-
- // Process variable.
-
//! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__process_variables}
auto const pv_config = config.getConfigSubtree("process_variables");
@@ -58,81 +55,48 @@ std::unique_ptr<Process> CreateTwoPhaseFlowWithPPProcess(
ProcessLib::parseSecondaryVariables(config, secondary_variables,
named_function_caller);
// Specific body force
- Eigen::VectorXd specific_body_force;
-
std::vector<double> const b =
//! \ogs_file_param{prj__processes__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());
+ Eigen::VectorXd specific_body_force(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{prj__processes__process__TWOPHASE_FLOW_PP__mass_lumping}
- auto mass_lump = config.getConfigParameter<bool>("mass_lumping");
+ auto const mass_lumping = config.getConfigParameter<bool>("mass_lumping");
+
+ auto& temperature = findParameter<double>(
+ config,
+ //! \ogs_file_param_special{prj__processes__process__TWOPHASE_FLOW_PP__temperature}
+ "temperature", parameters, 1);
//! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PP__material_property}
auto const& mat_config = config.getConfigSubtree("material_property");
- std::unique_ptr<
- MaterialLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPMaterialProperties>
- material = nullptr;
-
+ boost::optional<MeshLib::PropertyVector<int> const&> material_ids;
if (mesh.getProperties().existsPropertyVector<int>("MaterialIDs"))
{
INFO("The twophase flow is in heterogeneous porous media.");
- const bool has_material_ids = true;
- auto const& mat_ids =
- mesh.getProperties().getPropertyVector<int>("MaterialIDs");
- 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}};
+ material_ids =
+ *mesh.getProperties().getPropertyVector<int>("MaterialIDs");
}
else
{
INFO("The twophase flow is in homogeneous porous media.");
+ }
+ std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>
+ material = createTwoPhaseFlowWithPPMaterialProperties(mat_config, material_ids);
- MeshLib::Properties dummy_property;
-
- auto const& dummy_property_vector =
- dummy_property.createNewPropertyVector<int>(
- "MaterialIDs", MeshLib::MeshItemType::Cell, 1);
+ TwoPhaseFlowWithPPProcessData process_data{
+ specific_body_force, has_gravity, mass_lumping, temperature, std::move(material)};
- // 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}};
- }
+ 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
diff --git a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h
index 9f46501de39786dc09a40d1926e12442bba8d267..d6a43e27c3b2bb6b2ab6137b9ab2dc17f054fd06 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h
@@ -16,7 +16,7 @@ namespace ProcessLib
{
namespace TwoPhaseFlowWithPP
{
-std::unique_ptr<Process> CreateTwoPhaseFlowWithPPProcess(
+std::unique_ptr<Process> createTwoPhaseFlowWithPPProcess(
MeshLib::Mesh& mesh,
std::unique_ptr<ProcessLib::AbstractJacobianAssembler>&& jacobian_assembler,
std::vector<ProcessVariable> const& variables,
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
index 939b664468849ceb2f7f4bd12c507fcbbc59ec11..ee734e2f5e2f0a41acb0643c0e61f95d6c6b90c1 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler-impl.h
@@ -7,6 +7,27 @@
*
*/
+/**
+* common nomenclature
+* --------------primary variable----------------------
+* pn_int_pt pressure for nonwetting phase at each integration point
+* pc_int_pt capillary pressure at each integration point
+* --------------secondary variable--------------------
+* temperature capillary pressure
+* Sw wetting phase saturation
+* dSw_dpc derivative of wetting phase saturation with respect
+* to capillary pressure
+* rho_nonwet density of nonwetting phase
+* drhononwet_dpn derivative of nonwetting phase density with respect
+*to nonwetting phase pressure
+* rho_wet density of wetting phase
+* k_rel_nonwet relative permeability of nonwetting phase
+* mu_nonwet viscosity of nonwetting phase
+* lambda_nonwet mobility of nonwetting phase
+* k_rel_wet relative permeability of wetting phase
+* mu_wet viscosity of wetting phase
+* lambda_wet mobility of wetting phase
+*/
#pragma once
#include "TwoPhaseFlowWithPPLocalAssembler.h"
@@ -39,18 +60,12 @@ void TwoPhaseFlowWithPPLocalAssembler<
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);
@@ -61,9 +76,6 @@ void TwoPhaseFlowWithPPLocalAssembler<
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);
@@ -81,105 +93,97 @@ void TwoPhaseFlowWithPPLocalAssembler<
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)
+ const int material_id =
+ _process_data.material->getMaterialID(pos.getElementID().get());
+
+ const Eigen::MatrixXd& perm = _process_data.material->getPermeability(
+ material_id, t, pos, _element.getDimension());
+ assert(perm.rows() == _element.getDimension() || perm.rows() == 1);
+ GlobalDimMatrixType permeability = GlobalDimMatrixType::Zero(
+ _element.getDimension(), _element.getDimension());
+ if (perm.rows() == _element.getDimension())
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;
+ double pn_int_pt = 0.;
+ NumLib::shapeFunctionInterpolate(local_x, _ip_data[ip].N, pn_int_pt,
+ pc_int_pt);
- auto const& wp = _integration_method.getWeightedPoint(ip);
- const double integration_factor =
- sm.integralMeasure * sm.detJ * wp.getWeight();
+ _pressure_wet[ip] = pn_int_pt - pc_int_pt;
- double const rho_gas =
- _process_data._material->getGasDensity(pg_int_pt, _temperature);
- double const rho_w = _process_data._material->getLiquidDensity(
- _pressure_wetting[ip], _temperature);
+ const double temperature = _process_data.temperature(t, pos)[0];
+ double const rho_nonwet =
+ _process_data.material->getGasDensity(pn_int_pt, temperature);
+ double const rho_wet = _process_data.material->getLiquidDensity(
+ _pressure_wet[ip], temperature);
- double const Sw =
- (pc_int_pt < 0) ? 1.0 : interpolated_Pc.getValue(pc_int_pt);
+ double const Sw = _process_data.material->getSaturation(
+ material_id, t, pos, pn_int_pt, temperature, 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);
+ double dSw_dpc = _process_data.material->getSaturationDerivative(
+ material_id, t, pos, pn_int_pt, temperature, Sw);
+
+ double const porosity = _process_data.material->getPorosity(
+ material_id, t, pos, pn_int_pt, temperature, 0);
// Assemble M matrix
// nonwetting
- double const drhogas_dpg = _process_data._material->getDerivGasDensity(
- pg_int_pt, _temperature);
+ double const drhononwet_dpn =
+ _process_data.material->getGasDensityDerivative(pn_int_pt,
+ temperature);
- mass_operator.noalias() = sm.N.transpose() * sm.N * integration_factor;
+ Mgp.noalias() +=
+ porosity * (1 - Sw) * drhononwet_dpn * _ip_data[ip].massOperator;
+ Mgpc.noalias() +=
+ -porosity * rho_nonwet * dSw_dpc * _ip_data[ip].massOperator;
- 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;
+ Mlpc.noalias() +=
+ porosity * dSw_dpc * rho_wet * _ip_data[ip].massOperator;
- // 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;
+ double const k_rel_nonwet =
+ _process_data.material->getNonwetRelativePermeability(
+ t, pos, pn_int_pt, temperature, Sw);
+ double const mu_nonwet =
+ _process_data.material->getGasViscosity(pn_int_pt, temperature);
+ double const lambda_nonwet = k_rel_nonwet / mu_nonwet;
// 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;
+ double const k_rel_wet =
+ _process_data.material->getWetRelativePermeability(
+ t, pos, _pressure_wet[ip], temperature, Sw);
+ double const mu_wet = _process_data.material->getLiquidViscosity(
+ _pressure_wet[ip], temperature);
+ double const lambda_wet = k_rel_wet / mu_wet;
- laplace_operator.noalias() =
- sm.dNdx.transpose() * permeability * sm.dNdx * integration_factor;
+ laplace_operator.noalias() = _ip_data[ip].dNdx.transpose() *
+ permeability * _ip_data[ip].dNdx *
+ _ip_data[ip].integration_weight;
- 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;
+ Kgp.noalias() += rho_nonwet * lambda_nonwet * 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;
+ Klp.noalias() += rho_wet * lambda_wet * laplace_operator;
+ Klpc.noalias() += -rho_wet * lambda_wet * laplace_operator;
+ if (_process_data.has_gravity)
+ {
+ auto const& b = _process_data.specific_body_force;
+
+ NodalVectorType gravity_operator = _ip_data[ip].dNdx.transpose() *
+ permeability * b *
+ _ip_data[ip].integration_weight;
+ Bg.noalias() +=
+ rho_nonwet * rho_nonwet * lambda_nonwet * gravity_operator;
+ Bl.noalias() += rho_wet * rho_wet * lambda_wet * gravity_operator;
} // end of has gravity
- } // end of GP
- if (_process_data._has_mass_lumping)
+ }
+ if (_process_data.has_mass_lumping)
{
for (unsigned row = 0; row < Mgpc.cols(); row++)
{
@@ -193,8 +197,6 @@ void TwoPhaseFlowWithPPLocalAssembler<
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;
}
}
}
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
index f9d65553af82b05eaaee32a63fce4ec09ed7ad3f..2daf43b4cab64ed2e49b990a76f864160f3db467 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
@@ -20,13 +20,35 @@
#include "ProcessLib/Parameter/Parameter.h"
#include "ProcessLib/Utils/InitShapeMatrices.h"
-#include "MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h"
+#include "TwoPhaseFlowWithPPMaterialProperties.h"
#include "TwoPhaseFlowWithPPProcessData.h"
namespace ProcessLib
{
namespace TwoPhaseFlowWithPP
{
+template <typename NodalRowVectorType, typename GlobalDimNodalMatrixType,
+ typename NodalMatrixType>
+struct IntegrationPointData final
+{
+ explicit IntegrationPointData(
+ NodalRowVectorType const& N_, GlobalDimNodalMatrixType const& dNdx_,
+ TwoPhaseFlowWithPPMaterialProperties& material_property_,
+ double const& integration_weight_, NodalMatrixType const massOperator_)
+ : N(N_),
+ dNdx(dNdx_),
+ mat_property(material_property_),
+ integration_weight(integration_weight_),
+ massOperator(massOperator_)
+
+ {
+ }
+ NodalRowVectorType const N;
+ GlobalDimNodalMatrixType const dNdx;
+ TwoPhaseFlowWithPPMaterialProperties const& mat_property;
+ const double integration_weight;
+ NodalMatrixType const massOperator;
+};
const unsigned NUM_NODAL_DOF = 2;
class TwoPhaseFlowWithPPLocalAssemblerInterface
@@ -37,7 +59,7 @@ public:
virtual std::vector<double> const& getIntPtSaturation(
std::vector<double>& /*cache*/) const = 0;
- virtual std::vector<double> const& getIntPtWettingPressure(
+ virtual std::vector<double> const& getIntPtWetPressure(
std::vector<double>& /*cache*/) const = 0;
};
@@ -51,7 +73,10 @@ class TwoPhaseFlowWithPPLocalAssembler
using LocalAssemblerTraits = ProcessLib::LocalAssemblerTraits<
ShapeMatricesType, ShapeFunction::NPOINTS, NUM_NODAL_DOF, GlobalDim>;
+ using NodalRowVectorType = typename ShapeMatricesType::NodalRowVectorType;
+ using GlobalDimNodalMatrixType =
+ typename ShapeMatricesType::GlobalDimNodalMatrixType;
using NodalMatrixType = typename ShapeMatricesType::NodalMatrixType;
using NodalVectorType = typename ShapeMatricesType::NodalVectorType;
using GlobalDimMatrixType = typename ShapeMatricesType::GlobalDimMatrixType;
@@ -68,15 +93,29 @@ public:
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(
+ _pressure_wet(
std::vector<double>(_integration_method.getNumberOfPoints()))
{
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ _ip_data.reserve(n_integration_points);
+ auto const shape_matrices =
+ initShapeMatrices<ShapeFunction, ShapeMatricesType,
+ IntegrationMethod, GlobalDim>(
+ element, is_axially_symmetric, _integration_method);
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& sm = shape_matrices[ip];
+ _ip_data.emplace_back(
+ sm.N, sm.dNdx, *_process_data.material,
+ sm.integralMeasure * sm.detJ *
+ _integration_method.getWeightedPoint(ip).getWeight(),
+ sm.N.transpose() * sm.N * sm.integralMeasure * sm.detJ *
+ _integration_method.getWeightedPoint(ip).getWeight());
+ }
}
void assemble(double const t, std::vector<double> const& local_x,
@@ -87,7 +126,7 @@ public:
Eigen::Map<const Eigen::RowVectorXd> getShapeMatrix(
const unsigned integration_point) const override
{
- auto const& N = _shape_matrices[integration_point].N;
+ auto const& N = _ip_data[integration_point].N;
// assumes N is stored contiguously in memory
return Eigen::Map<const Eigen::RowVectorXd>(N.data(), N.size());
@@ -100,28 +139,32 @@ public:
return _saturation;
}
- std::vector<double> const& getIntPtWettingPressure(
+ std::vector<double> const& getIntPtWetPressure(
std::vector<double>& /*cache*/) const override
{
- assert(_pressure_wetting.size() > 0);
- return _pressure_wetting;
+ assert(_pressure_wet.size() > 0);
+ return _pressure_wet;
}
private:
MeshLib::Element const& _element;
IntegrationMethod const _integration_method;
- std::vector<ShapeMatrices, Eigen::aligned_allocator<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<
+ IntegrationPointData<NodalRowVectorType, GlobalDimNodalMatrixType,
+ NodalMatrixType>,
+ Eigen::aligned_allocator<IntegrationPointData<
+ NodalRowVectorType, GlobalDimNodalMatrixType, NodalMatrixType>>>
+ _ip_data;
+
+ // output vector for wetting phase saturation with
+ // respect to each integration point
std::vector<double> _saturation;
- std::vector<double> _pressure_wetting;
+ // output vector for wetting phase pressure with respect
+ // to each integration point
+ std::vector<double> _pressure_wet;
static const int nonwet_pressure_coeff_index = 0;
static const int cap_pressure_coeff_index = 1;
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c2ed90e5a791f7822116615a74ad4eaf2df7e688
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.cpp
@@ -0,0 +1,169 @@
+/**
+* \copyright
+* Copyright (c) 2012-2017, 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 "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.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 ProcessLib
+{
+namespace TwoPhaseFlowWithPP
+{
+TwoPhaseFlowWithPPMaterialProperties::TwoPhaseFlowWithPPMaterialProperties(
+ boost::optional<MeshLib::PropertyVector<int> const&> const material_ids,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ liquid_density,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ liquid_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,
+ std::vector<std::unique_ptr<
+ MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
+ capillary_pressure_models,
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
+ relative_permeability_models)
+ : _liquid_density(std::move(liquid_density)),
+ _liquid_viscosity(std::move(liquid_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)),
+ _capillary_pressure_models(std::move(capillary_pressure_models)),
+ _relative_permeability_models(std::move(relative_permeability_models))
+{
+ DBUG("Create material properties for Two-Phase flow with PP model.");
+}
+
+int TwoPhaseFlowWithPPMaterialProperties::getMaterialID(
+ const std::size_t element_id)
+{
+ if (!_material_ids)
+ {
+ return 0;
+ }
+
+ assert(element_id < _material_ids->size());
+ return (*_material_ids)[element_id];
+}
+
+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::p)] = 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::p)] = p;
+ return _gas_density->getValue(vars);
+}
+
+double TwoPhaseFlowWithPPMaterialProperties::getGasDensityDerivative(
+ 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::p)] = p;
+
+ return _gas_density->getdValue(vars,
+ MaterialLib::Fluid::PropertyVariableType::p);
+}
+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::p)] = p;
+ return _liquid_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::p)] = p;
+ return _gas_viscosity->getValue(vars);
+}
+
+Eigen::MatrixXd const& TwoPhaseFlowWithPPMaterialProperties::getPermeability(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const int /*dim*/) const
+{
+ return _intrinsic_permeability_models[material_id];
+}
+
+double TwoPhaseFlowWithPPMaterialProperties::getPorosity(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+ const double T, const double porosity_variable) const
+{
+ return _porosity_models[material_id]->getValue(porosity_variable, T);
+}
+
+double TwoPhaseFlowWithPPMaterialProperties::getNonwetRelativePermeability(
+ const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/,
+ const double /*p*/, const double /*T*/, const double saturation) const
+{
+ return _relative_permeability_models[0]->getValue(saturation);
+}
+
+double TwoPhaseFlowWithPPMaterialProperties::getWetRelativePermeability(
+ const double /*t*/, const ProcessLib::SpatialPosition& /*pos*/,
+ const double /*p*/, const double /*T*/, const double saturation) const
+{
+ return _relative_permeability_models[1]->getValue(saturation);
+}
+
+double TwoPhaseFlowWithPPMaterialProperties::getSaturation(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+ const double /*T*/, const double pc) const
+{
+ return _capillary_pressure_models[material_id]->getSaturation(pc);
+}
+double TwoPhaseFlowWithPPMaterialProperties::getSaturationDerivative(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+ const double /*T*/, const double saturation) const
+{
+ const double dpcdsw =
+ _capillary_pressure_models[material_id]->getdPcdS(saturation);
+ return 1 / dpcdsw;
+}
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..fd3b0657d019e9928c962bd6474bcb2743c4f11d
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPMaterialProperties.h
@@ -0,0 +1,129 @@
+/**
+* \copyright
+* Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+* Distributed under a Modified BSD License.
+* See accompanying file LICENSE.txt or
+* http://www.opengeosys.org/project/license
+*
+*/
+
+#pragma once
+
+#include <memory>
+#include <vector>
+#include "MaterialLib/Fluid/FluidPropertyHeaders.h"
+#include "MaterialLib/PhysicalConstant.h"
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/PorousPropertyHeaders.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/CreateRelativePermeabilityModel.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/RelativePermeability/RelativePermeability.h"
+
+namespace MeshLib
+{
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
+namespace ProcessLib
+{
+class SpatialPosition;
+namespace TwoPhaseFlowWithPP
+{
+ /** This class has a collection of material properties for two-phase flow with PP model
+ * and it makes description of the material properties for two-phase condition,
+ * i.e. the gas/liquid density and viscosity models, respectively,
+ * the relative permeability models with respect to two phases,
+ * the capillary pressure-saturation relationships.
+ * It generally provides the computation of the PDE coefficients for two-phase flow.
+ */
+
+class TwoPhaseFlowWithPPMaterialProperties
+{
+public:
+ using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType;
+
+ TwoPhaseFlowWithPPMaterialProperties(
+ boost::optional<MeshLib::PropertyVector<int> const&> const material_ids,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ liquid_density,
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty>
+ liquid_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,
+ std::vector<std::unique_ptr<
+ MaterialLib::PorousMedium::CapillaryPressureSaturation>>&&
+ capillary_pressure_models,
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>&&
+ relative_permeability_models);
+
+ int getMaterialID(const std::size_t element_id);
+
+ Eigen::MatrixXd const& getPermeability(
+ const int material_id,
+ const double t,
+ const ProcessLib::SpatialPosition& pos,
+ const int dim) const;
+
+ double getPorosity(const int material_id, const double t,
+ const ProcessLib::SpatialPosition& pos, const double p,
+ const double T, const double porosity_variable) const;
+
+ double getNonwetRelativePermeability(const double t,
+ const ProcessLib::SpatialPosition& pos,
+ const double p, const double T,
+ const double saturation) const;
+ double getWetRelativePermeability(const double t,
+ const ProcessLib::SpatialPosition& pos,
+ const double p, const double T,
+ const double saturation) const;
+ double getSaturation(const int material_id, const double t,
+ const ProcessLib::SpatialPosition& pos, const double p,
+ const double T, const double pc) const;
+ double getSaturationDerivative(const int material_id, const double t,
+ const ProcessLib::SpatialPosition& pos,
+ const double p, const double T,
+ const double saturation) 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 getGasDensityDerivative(double const p, double const T) const;
+
+protected:
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_density;
+ std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_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.
+ */
+ boost::optional<MeshLib::PropertyVector<int> const&> const _material_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;
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::CapillaryPressureSaturation>>
+ _capillary_pressure_models;
+ std::vector<
+ std::unique_ptr<MaterialLib::PorousMedium::RelativePermeability>>
+ _relative_permeability_models;
+};
+
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
index 2cbb77f5b51e0b3f963aa80508374ba02e887c31..6c855e3ba4cb4ce6e64cea3ecbdfb6ad9e519736 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.cpp
@@ -35,7 +35,7 @@ TwoPhaseFlowWithPPProcess::TwoPhaseFlowWithPPProcess(
BaseLib::ConfigTree const& /*config*/,
std::map<std::string,
std::unique_ptr<MathLib::PiecewiseLinearInterpolation>> const&
- /*curves*/)
+ /*curves*/)
: Process(mesh, std::move(jacobian_assembler), parameters,
integration_order, std::move(process_variables),
std::move(secondary_variables), std::move(named_function_caller)),
@@ -62,10 +62,10 @@ void TwoPhaseFlowWithPPProcess::initializeConcreteProcess(
&TwoPhaseFlowWithPPLocalAssemblerInterface::getIntPtSaturation));
_secondary_variables.addSecondaryVariable(
- "pressure_wetting", 1,
+ "pressure_wet", 1,
makeExtrapolator(getExtrapolator(), _local_assemblers,
&TwoPhaseFlowWithPPLocalAssemblerInterface::
- getIntPtWettingPressure));
+ getIntPtWetPressure));
}
void TwoPhaseFlowWithPPProcess::assembleConcreteProcess(const double t,
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
index 42854353fe426ccf8c388066ffe5a8d900cec95f..376d91d356f1ed5823af425e1234b861c1ae065a 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcess.h
@@ -9,7 +9,6 @@
#pragma once
-#include "MaterialLib/TwoPhaseModels/TwoPhaseFlowWithPPMaterialProperties.h"
#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
#include "NumLib/DOF/LocalToGlobalIndexMap.h"
#include "ProcessLib/Process.h"
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h
index 6da20a21694f2252fd2427e17b290ceef3a10a1b..bf722b2473d66603b603ddd794f73c62a29181e0 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPProcessData.h
@@ -8,7 +8,7 @@
*/
#pragma once
-
+#include "TwoPhaseFlowWithPPMaterialProperties.h"
namespace MeshLib
{
class Element;
@@ -27,29 +27,23 @@ struct 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_)
+ Parameter<double> const& temperature_,
+ std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties>&& material_)
+ : specific_body_force(specific_body_force_),
+ has_gravity(has_gravity_),
+ has_mass_lumping(has_mass_lumping_),
+ temperature(temperature_),
+ material(std::move(material_))
+
{
}
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)
+ : specific_body_force(other.specific_body_force),
+ has_gravity(other.has_gravity),
+ has_mass_lumping(other.has_mass_lumping),
+ temperature(other.temperature),
+ material(std::move(other.material))
{
}
@@ -66,18 +60,14 @@ struct TwoPhaseFlowWithPPProcessData
//! Specific body forces applied to solid and fluid.
//! It is usually used to apply gravitational forces.
//! A vector of displacement dimension's length.
- Eigen::VectorXd const _specific_body_force;
+ Eigen::VectorXd const specific_body_force;
- bool const _has_gravity;
+ bool const has_gravity;
//! Enables lumping of the mass matrix.
- 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;
+ bool const has_mass_lumping;
+ Parameter<double> const& temperature;
+ std::unique_ptr<TwoPhaseFlowWithPPMaterialProperties> material;
};
} // namespace TwoPhaseFlowWithPP