diff --git a/Applications/ApplicationsLib/ProjectData.cpp b/Applications/ApplicationsLib/ProjectData.cpp
index 08be9bd3f86526061e50163646f42bceeee920ce..2436b1a083b93c709cd2d1e725ae511e0c878735 100644
--- a/Applications/ApplicationsLib/ProjectData.cpp
+++ b/Applications/ApplicationsLib/ProjectData.cpp
@@ -35,16 +35,17 @@
#include "ProcessLib/UncoupledProcessesTimeLoop.h"
#include "ProcessLib/GroundwaterFlow/CreateGroundwaterFlowProcess.h"
+#include "ProcessLib/HT/CreateHTProcess.h"
#include "ProcessLib/HeatConduction/CreateHeatConductionProcess.h"
#include "ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.h"
#include "ProcessLib/LIE/HydroMechanics/CreateHydroMechanicsProcess.h"
#include "ProcessLib/LIE/SmallDeformation/CreateSmallDeformationProcess.h"
-#include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
#include "ProcessLib/LiquidFlow/CreateLiquidFlowProcess.h"
+#include "ProcessLib/RichardsFlow/CreateRichardsFlowProcess.h"
#include "ProcessLib/SmallDeformation/CreateSmallDeformationProcess.h"
#include "ProcessLib/TES/CreateTESProcess.h"
-#include "ProcessLib/HT/CreateHTProcess.h"
#include "ProcessLib/TwoPhaseFlowWithPP/CreateTwoPhaseFlowWithPPProcess.h"
+#include "ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowWithPrhoProcess.h"
namespace detail
{
@@ -442,6 +443,15 @@ void ProjectData::parseProcesses(BaseLib::ConfigTree const& processes_config,
process_config, _curves);
}
+ else if (type == "TWOPHASE_FLOW_PRHO")
+ {
+ process = ProcessLib::TwoPhaseFlowWithPrho::
+ createTwoPhaseFlowWithPrhoProcess(
+ *_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());
@@ -526,8 +536,8 @@ void ProjectData::parseCurves(
BaseLib::insertIfKeyUniqueElseError(
_curves,
name,
- MathLib::createPiecewiseLinearCurve
- <MathLib::PiecewiseLinearInterpolation>(conf),
+ MathLib::createPiecewiseLinearCurve<
+ MathLib::PiecewiseLinearInterpolation>(conf),
"The curve name is not unique.");
}
}
diff --git a/Documentation/ProjectFile/material/porous_medium/capillary_pressure/BrookCorey/t_sg_r.md b/Documentation/ProjectFile/material/porous_medium/capillary_pressure/BrookCorey/t_sg_r.md
new file mode 100644
index 0000000000000000000000000000000000000000..cbc92bb673eb1076b8590c712a5fefcc6b5fc892
--- /dev/null
+++ b/Documentation/ProjectFile/material/porous_medium/capillary_pressure/BrookCorey/t_sg_r.md
@@ -0,0 +1 @@
+A tag for the residual saturation in the gas phase.
diff --git a/Documentation/ProjectFile/material/porous_medium/capillary_pressure/vanGenuchten/t_has_regularized.md b/Documentation/ProjectFile/material/porous_medium/capillary_pressure/vanGenuchten/t_has_regularized.md
new file mode 100644
index 0000000000000000000000000000000000000000..9627930e3d621f8098d03e6c446f7f29ecc3536d
--- /dev/null
+++ b/Documentation/ProjectFile/material/porous_medium/capillary_pressure/vanGenuchten/t_has_regularized.md
@@ -0,0 +1 @@
+A tag for using regularized van Genuchten model.
diff --git a/Documentation/ProjectFile/material/porous_medium/capillary_pressure/vanGenuchten/t_sg_r.md b/Documentation/ProjectFile/material/porous_medium/capillary_pressure/vanGenuchten/t_sg_r.md
new file mode 100644
index 0000000000000000000000000000000000000000..b8e8041fe4c8cc63dbf4d704847ca4ce6640cd58
--- /dev/null
+++ b/Documentation/ProjectFile/material/porous_medium/capillary_pressure/vanGenuchten/t_sg_r.md
@@ -0,0 +1 @@
+A tag for the residual saturation of the nonwetting phase.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/c_TWOPHASE_FLOW_PRHO.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/c_TWOPHASE_FLOW_PRHO.md
new file mode 100644
index 0000000000000000000000000000000000000000..af0e06cd87dec000f83dd52326b6a9fe3f888896
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/c_TWOPHASE_FLOW_PRHO.md
@@ -0,0 +1 @@
+A tag for compositional based two-phase flow process, using liquid pressure and total mass density as process variables.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/i_material_property.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/i_material_property.md
new file mode 100644
index 0000000000000000000000000000000000000000..af3a673dfd6b59c235f84676a32ddc6fed46564b
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/i_material_property.md
@@ -0,0 +1 @@
+A tag for liquid and gas material properties of the two-phase flow process.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/i_porous_medium.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/i_porous_medium.md
new file mode 100644
index 0000000000000000000000000000000000000000..dd17a8239ab81825b29d17d3562b192a40286e28
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/i_porous_medium.md
@@ -0,0 +1 @@
+A tag for the properties of porous media.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/a_id.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/a_id.md
new file mode 100644
index 0000000000000000000000000000000000000000..dcddf575a6e74058d49ad1b156ae63d21829778c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/a_id.md
@@ -0,0 +1 @@
+It specifies the material ID of a porous medium. The material IDs of elements should be given in the mesh's cell data array "MaterialIDs".
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/i_porous_medium.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/i_porous_medium.md
new file mode 100644
index 0000000000000000000000000000000000000000..eb12008f2e28c580da129187cd9f843209a75930
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/i_porous_medium.md
@@ -0,0 +1 @@
+A tag for the properties of a porous medium with material ID.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/relative_permeability/i_relative_permeability.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/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_PRHO/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_PRHO/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/a_id.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/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_PRHO/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_PRHO/material_property/porous_medium/porous_medium/relative_permeability/relative_permeability/i_relative_permeability.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/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_PRHO/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_PRHO/material_property/porous_medium/porous_medium/t_capillary_pressure.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/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_PRHO/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_PRHO/material_property/porous_medium/porous_medium/t_permeability.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_permeability.md
new file mode 100644
index 0000000000000000000000000000000000000000..0f0e1d4fb6be8191211eb35bc6d5abee1065e5ff
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_permeability.md
@@ -0,0 +1 @@
+A tag for the intrinsic permeability model.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_porosity.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_porosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..45e8e708f526fc6bd69b42ac3cc5e19ee0a07d8c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_porosity.md
@@ -0,0 +1 @@
+A tag for the porosity model.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_storage.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_storage.md
new file mode 100644
index 0000000000000000000000000000000000000000..f036cada72324684a9a256be78cfb647ed8b5fe6
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/porous_medium/porous_medium/t_storage.md
@@ -0,0 +1 @@
+A tag for the storage model.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_fluid.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_fluid.md
new file mode 100644
index 0000000000000000000000000000000000000000..6f73cf4d6bfb9b1b683522a9a996b689e594cd06
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_fluid.md
@@ -0,0 +1 @@
+Defines fluid properties of gas and liquid phases.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_gas_density.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_gas_density.md
new file mode 100644
index 0000000000000000000000000000000000000000..bcf97b55a9e2bf33941d089ec03b97f90d562b1f
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/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_PRHO/material_property/t_gas_viscosity.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_gas_viscosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..a02e98700f601132ce1e4e545b1f87d2d0eb5388
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/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_PRHO/material_property/t_liquid_density.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_liquid_density.md
new file mode 100644
index 0000000000000000000000000000000000000000..5c8db5dec3d6f9cc8f530470cbca94b43e23411c
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/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_PRHO/material_property/t_liquid_viscosity.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/material_property/t_liquid_viscosity.md
new file mode 100644
index 0000000000000000000000000000000000000000..0451570185b81e26ecab7dad8eb9b7af93776891
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/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_PRHO/process_variables/i_process_variables.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/i_process_variables.md
new file mode 100644
index 0000000000000000000000000000000000000000..530b3234e30086980d5c5b39edb24619316c0552
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/i_process_variables.md
@@ -0,0 +1 @@
+The process variables for liquid pressure and overall mass density.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/t_liquid_pressure.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/t_liquid_pressure.md
new file mode 100644
index 0000000000000000000000000000000000000000..fd5abc31467988953fe4b78734ba0b17af757ecf
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/t_liquid_pressure.md
@@ -0,0 +1 @@
+Pressure for the liquid phase.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/t_overall_mass_density.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/t_overall_mass_density.md
new file mode 100644
index 0000000000000000000000000000000000000000..43608935aaad5648523240c77ee9b2665dcba08a
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/process_variables/t_overall_mass_density.md
@@ -0,0 +1 @@
+Overall mass density of the light component
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_diffusion_coeff_component_a.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_diffusion_coeff_component_a.md
new file mode 100644
index 0000000000000000000000000000000000000000..dcc47ec43e37af1e141c753359fdf9a23c606796
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_diffusion_coeff_component_a.md
@@ -0,0 +1 @@
+Diffusion coefficient for the heavy component.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_diffusion_coeff_component_b.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_diffusion_coeff_component_b.md
new file mode 100644
index 0000000000000000000000000000000000000000..2953459bc7b1a1fbb446b0cacca6fe969c973a10
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_diffusion_coeff_component_b.md
@@ -0,0 +1 @@
+Diffusion coefficient for the light component.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_mass_lumping.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_mass_lumping.md
new file mode 100644
index 0000000000000000000000000000000000000000..51c8c6682bb70b7ccb94aa722519fa4b9c9b2363
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_mass_lumping.md
@@ -0,0 +1 @@
+A tag for using mass-lumping technique.
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_specific_body_force.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_specific_body_force.md
new file mode 100644
index 0000000000000000000000000000000000000000..9adb4107b6da4bb03a2d35254e56812414eeb74f
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_specific_body_force.md
@@ -0,0 +1 @@
+\copydoc ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoProcessData::_specific_body_force
diff --git a/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_temperature.md b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_temperature.md
new file mode 100644
index 0000000000000000000000000000000000000000..bf468ba780d421455d62231de19de8c75a1a1211
--- /dev/null
+++ b/Documentation/ProjectFile/prj/processes/process/TWOPHASE_FLOW_PRHO/t_temperature.md
@@ -0,0 +1 @@
+Constant temperature (K) for the simulation.
diff --git a/MaterialLib/PhysicalConstant.h b/MaterialLib/PhysicalConstant.h
index cce50a0e92d9a65b234a8ff4038ce0eec58f7c9d..47ddab80ee36c62c8418ac4a1044c445552cf0de 100644
--- a/MaterialLib/PhysicalConstant.h
+++ b/MaterialLib/PhysicalConstant.h
@@ -74,6 +74,13 @@ const double O2 = 0.032; ///< kg \per{mol}
* - Source: http://www.engineeringtoolbox.com/molecular-mass-air-d_679.html
*/
const double Air = 0.02897; ///< kg \per{mol}
+
+/**
+ * Hydrogen.
+ *
+ * - Source: https://pubchem.ncbi.nlm.nih.gov/compound/Hydrogen
+ */
+const double H2 = 0.002016; ///< kg \per{mol}
} // namespace MolarMass
/**
@@ -86,5 +93,18 @@ namespace SpecificGasConstant
/// Specific gas constant for water vapour.
const double WaterVapour = IdealGasConstant / MolarMass::Water; // = 461.504;
} // namespace SpecificGasConstant
+/**
+* Henry's law constant
+* It is defined here as the ratio of the aqueous-phase concentration of a
+* chemical to its equilibrium partial pressure in the gas phase.
+*/
+namespace HenryConstant
+{
+/**Henry constant for hydrogen, Please refer to
+* De Nevers N. Physical and chemical equilibrium for chemical engineers[M].
+* John Wiley & Sons, 2012.
+*/
+double const HenryConstantH2 = 7.65e-6; /// mol/Pa./m3
+} // namespace HenryConstant
} // namespace PhysicalConstant
} // namespace MaterialLib
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h
index 519f2ee1a9687acb4733731b46599ebcbcf22ac4..f07b5af7ba0503f3742b12f7ec50eb88dd70b50a 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h
@@ -44,14 +44,15 @@ public:
/**
* @param pb Entry pressure, \f$ p_b \f$
* @param Sr Residual saturation, \f$ S_r \f$
+ * @param Sg_r Residual saturation of gas phase, \f$ Sg_r \f$
* @param Smax Maximum saturation, \f$ S_{\mbox{max}} \f$
* @param m Exponent, \f$ m \f$
* @param Pc_max Maximum capillary pressure, \f$ P_c^{\mbox{max}}\f$
*/
BrookCoreyCapillaryPressureSaturation(const double pb, const double Sr,
- const double Smax, const double m,
- const double Pc_max)
- : CapillaryPressureSaturation(Sr, Smax, Pc_max), _pb(pb), _m(m)
+ const double Sg_r, const double Smax,
+ const double m, const double Pc_max)
+ : CapillaryPressureSaturation(Sr, Sg_r, Smax, Pc_max), _pb(pb), _m(m)
{
}
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h
index b3dbf9dc6ff58625c19b777316b30ddade7971e9..04a6b7b1175d5c89a878b2b3cfa67581e70ebf8a 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h
@@ -23,12 +23,16 @@ class CapillaryPressureSaturation
public:
/**
* @param Sr Residual saturation, \f$ S_r \f$
+ * @param Sg_r Residual saturation of gas phase, \f$ Sg_r \f$
* @param Smax Maximum saturation, \f$ S_{\mbox{max}} \f$
* @param Pc_max Maximum capillary pressure, \f$ P_c^{\mbox{max}}\f$
*/
- CapillaryPressureSaturation(const double Sr, const double Smax,
- const double Pc_max)
- : _saturation_r(Sr), _saturation_max(Smax), _pc_max(Pc_max)
+ CapillaryPressureSaturation(const double Sr, const double Sg_r,
+ const double Smax, const double Pc_max)
+ : _saturation_r(Sr),
+ _saturation_nonwet_r(Sg_r),
+ _saturation_max(Smax),
+ _pc_max(Pc_max)
{
}
virtual ~CapillaryPressureSaturation() = default;
@@ -38,7 +42,6 @@ public:
/// Get capillary pressure.
virtual double getCapillaryPressure(const double saturation) const = 0;
-
/// Get saturation.
virtual double getSaturation(const double capillary_ressure) const = 0;
@@ -46,9 +49,11 @@ public:
virtual double getdPcdS(const double saturation) const = 0;
protected:
- const double _saturation_r; ///< Residual saturation.
- const double _saturation_max; ///< Maximum saturation.
- const double _pc_max; ///< Maximum capillaray pressure
+ const double _saturation_r; ///< Residual saturation.
+ const double _saturation_nonwet_r; ///< Residual saturation of nonwetting
+ ///phase (optional).
+ const double _saturation_max; ///< Maximum saturation.
+ const double _pc_max; ///< Maximum capillaray pressure
/** A small number for an offset:
* 1. to set the bound of S, the saturation, such that
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturationCurve.h b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturationCurve.h
index 63dc649c1956cf787a4f6647ca35d780e21078c4..fd55b5ac0d547cf82f374665d44a07110f8e3388 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturationCurve.h
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturationCurve.h
@@ -31,6 +31,7 @@ public:
std::unique_ptr<MathLib::PiecewiseLinearMonotonicCurve>&& curve_data)
: CapillaryPressureSaturation(
curve_data->getSupportMin(),
+ 1.0 - curve_data->getSupportMax(),
curve_data->getSupportMax(),
curve_data->getValue(curve_data->getSupportMin())),
_curve_data(std::move(curve_data))
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp
index 7a25f0fadfb29d8574132ef01235414b1cf09c19..ecbb34315b1ef76a25b44800e2e378f19fb4d546 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp
@@ -44,6 +44,16 @@ static std::unique_ptr<CapillaryPressureSaturation> createBrookCorey(
//! \ogs_file_param{material__porous_medium__capillary_pressure__BrookCorey__sr}
const double Sr = config.getConfigParameter<double>("sr");
+ double Sg_r = 0.0;
+ //! \ogs_file_param{material__porous_medium__capillary_pressure__BrookCorey__sg_r}
+ if (auto const Sg_r_ptr = config.getConfigParameterOptional<double>("sg_r"))
+ {
+ DBUG(
+ "Using value %g for nonwetting phase residual saturation in "
+ "capillary pressure model.",
+ (*Sg_r_ptr));
+ Sg_r = *Sg_r_ptr;
+ }
//! \ogs_file_param{material__porous_medium__capillary_pressure__BrookCorey__smax}
const double Smax = config.getConfigParameter<double>("smax");
@@ -59,7 +69,8 @@ static std::unique_ptr<CapillaryPressureSaturation> createBrookCorey(
const double Pc_max = config.getConfigParameter<double>("pc_max");
return std::unique_ptr<CapillaryPressureSaturation>(
- new BrookCoreyCapillaryPressureSaturation(pd, Sr, Smax, m, Pc_max));
+ new BrookCoreyCapillaryPressureSaturation(
+ pd, Sr, Sg_r, Smax, m, Pc_max));
}
/**
@@ -79,6 +90,17 @@ static std::unique_ptr<CapillaryPressureSaturation> createVanGenuchten(
//! \ogs_file_param{material__porous_medium__capillary_pressure__vanGenuchten__sr}
const double Sr = config.getConfigParameter<double>("sr");
+ double Sg_r = 0.0;
+ //! \ogs_file_param{material__porous_medium__capillary_pressure__vanGenuchten__sg_r}
+ if (auto const Sg_r_ptr = config.getConfigParameterOptional<double>("sg_r"))
+ {
+ DBUG(
+ "Using value %g for nonwetting phase residual saturation in "
+ "capillary pressure model.",
+ (*Sg_r_ptr));
+ Sg_r = *Sg_r_ptr;
+ }
+
//! \ogs_file_param{material__porous_medium__capillary_pressure__vanGenuchten__smax}
const double Smax = config.getConfigParameter<double>("smax");
@@ -93,8 +115,18 @@ static std::unique_ptr<CapillaryPressureSaturation> createVanGenuchten(
//! \ogs_file_param{material__porous_medium__capillary_pressure__vanGenuchten__pc_max}
const double Pc_max = config.getConfigParameter<double>("pc_max");
+ bool has_regularized = false;
+ if (auto const has_regularized_conf =
+ //! \ogs_file_param{material__porous_medium__capillary_pressure__vanGenuchten__has_regularized}
+ config.getConfigParameterOptional<bool>("has_regularized"))
+ {
+ DBUG("capillary pressure model: %s",
+ (*has_regularized_conf) ? "true" : "false");
+ has_regularized = *has_regularized_conf;
+ }
return std::unique_ptr<CapillaryPressureSaturation>(
- new VanGenuchtenCapillaryPressureSaturation(pd, Sr, Smax, m, Pc_max));
+ new VanGenuchtenCapillaryPressureSaturation(
+ pd, Sr, Sg_r, Smax, m, Pc_max, has_regularized));
}
std::unique_ptr<CapillaryPressureSaturation> createCapillaryPressureModel(
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.cpp
index 924ae326e26684b079cb85529d468ac693df0189..5e7e028a7cad49575ab4fb189252062fe20ed97a 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.cpp
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.cpp
@@ -23,6 +23,25 @@ namespace PorousMedium
double VanGenuchtenCapillaryPressureSaturation::getCapillaryPressure(
const double saturation) const
{
+ if (_has_regularized)
+ {
+ double Sg = 1 - saturation;
+ if (Sg <= 1 - _saturation_r && Sg >= _saturation_nonwet_r)
+ {
+ return getPcBarvGSg(Sg);
+ }
+ else if (Sg < _saturation_nonwet_r)
+ {
+ return getPcBarvGSg(_saturation_nonwet_r) +
+ getdPcdSvGBar(_saturation_nonwet_r) *
+ (Sg - _saturation_nonwet_r);
+ }
+ else
+ {
+ return getPcBarvGSg(1 - _saturation_r) +
+ getdPcdSvGBar(1 - _saturation_r) * (Sg - 1 + _saturation_r);
+ }
+ }
const double S =
MathLib::limitValueInInterval(saturation, _saturation_r + _minor_offset,
_saturation_max - _minor_offset);
@@ -46,6 +65,22 @@ double VanGenuchtenCapillaryPressureSaturation::getSaturation(
double VanGenuchtenCapillaryPressureSaturation::getdPcdS(
const double saturation) const
{
+ if (_has_regularized)
+ {
+ double const Sg = 1 - saturation;
+ if (Sg >= _saturation_nonwet_r && Sg <= 1 - _saturation_r)
+ {
+ return -getdPcdSvGBar(Sg);
+ }
+ else if (Sg < _saturation_nonwet_r)
+ {
+ return -getdPcdSvGBar(_saturation_nonwet_r);
+ }
+ else
+ {
+ return -getdPcdSvGBar(1 - _saturation_r);
+ }
+ }
const double S =
MathLib::limitValueInInterval(saturation, _saturation_r + _minor_offset,
_saturation_max - _minor_offset);
@@ -54,6 +89,50 @@ double VanGenuchtenCapillaryPressureSaturation::getdPcdS(
const double val2 = std::pow(val1 - 1.0, -_m);
return _pb * (_m - 1.0) * val1 * val2 / (_m * (S - _saturation_r));
}
+/// Regularized van Genuchten capillary pressure-saturation Model
+double VanGenuchtenCapillaryPressureSaturation::getPcBarvGSg(double Sg) const
+{
+ double const Sg_r = CapillaryPressureSaturation::_saturation_nonwet_r;
+ double const S_lr = CapillaryPressureSaturation::_saturation_r;
+ double const S_bar = getSBar(Sg);
+ return getPcvGSg(S_bar) - getPcvGSg(Sg_r + (1 - Sg_r - S_lr) * _xi / 2);
+}
+/// Regularized van Genuchten capillary pressure-saturation Model
+double VanGenuchtenCapillaryPressureSaturation::getSBar(double Sg) const
+{
+ double const Sg_r = CapillaryPressureSaturation::_saturation_nonwet_r;
+ double const S_lr = CapillaryPressureSaturation::_saturation_r;
+ return Sg_r + (1 - _xi) * (Sg - Sg_r) + 0.5 * _xi * (1 - Sg_r - S_lr);
+}
+/// van Genuchten capillary pressure-saturation Model
+double VanGenuchtenCapillaryPressureSaturation::getPcvGSg(double Sg) const
+{
+ double const Sg_r = CapillaryPressureSaturation::_saturation_nonwet_r;
+ double const S_lr = CapillaryPressureSaturation::_saturation_r;
+ double const S_le = (1 - Sg - S_lr) /
+ (1 - Sg_r - CapillaryPressureSaturation::_saturation_r);
+ return _pb * std::pow(std::pow(S_le, (-1.0 / _m)) - 1.0, 1.0 - _m);
+}
+/// derivative dPCdS based on regularized van Genuchten capillary
+/// pressure-saturation Model
+double VanGenuchtenCapillaryPressureSaturation::getdPcdSvGBar(double Sg) const
+{
+ double S_bar = getSBar(Sg);
+ return getdPcdSvG(S_bar) * (1 - _xi);
+}
+/// derivative dPCdS based on standard van Genuchten capillary
+/// pressure-saturation Model
+double VanGenuchtenCapillaryPressureSaturation::getdPcdSvG(
+ const double Sg) const
+{
+ double const Sg_r = CapillaryPressureSaturation::_saturation_nonwet_r;
+ double const S_lr = CapillaryPressureSaturation::_saturation_r;
+ double const nn = 1 / (1 - _m);
+ double const S_le = (1 - Sg - S_lr) / (1 - Sg_r - S_lr);
+ return _pb * (1 / (_m * nn)) * (1 / (1 - S_lr - Sg_r)) *
+ std::pow(std::pow(S_le, (-1 / _m)) - 1, (1 / nn) - 1) *
+ std::pow(S_le, (-1 / _m)) / S_le;
+}
} // end namespace
} // end namespace
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.h b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.h
index 3c60b0d1fa42a04c0b8ba2a8824c18870afe5e27..b791e1334b6eec3a71cb37f43c53a893243ffe5a 100644
--- a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.h
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.h
@@ -12,8 +12,8 @@
#pragma once
+#include <cmath>
#include <limits>
-
#include "CapillaryPressureSaturation.h"
namespace MaterialLib
@@ -42,6 +42,14 @@ namespace PorousMedium
* as
* \f[p_b=\rho g/\alpha\f]
*/
+/**
+* \brief The regularized van Genuchten model, please ref to
+* Marchand E, Mueller T, Knabner P.
+* Fully coupled generalized hybrid-mixed finite element approximation of
+* two-phase two-component flow in porous media.
+* Part I: formulation and properties of the mathematical model. Comput Geosci
+* 2013;17(2):431每42.
+*/
class VanGenuchtenCapillaryPressureSaturation final
: public CapillaryPressureSaturation
{
@@ -49,14 +57,22 @@ public:
/**
* @param pb Entry pressure, \f$ p_b \f$
* @param Sr Residual saturation, \f$ S_r \f$
+ * @param Sg_r Residual saturation, \f$ S_g^{\mbox{r}} \f$
* @param Smax Maximum saturation, \f$ S_{\mbox{max}} \f$
* @param m Exponent, \f$ m \f$
* @param Pc_max Maximum capillary pressure, \f$ P_c^{\mbox{max}}\f$
+ * @param has_regularized whether use the regularized van Genuchten model,
+ * \f$ m \f$
*/
VanGenuchtenCapillaryPressureSaturation(const double pb, const double Sr,
+ const double Sg_r,
const double Smax, const double m,
- const double Pc_max)
- : CapillaryPressureSaturation(Sr, Smax, Pc_max), _pb(pb), _m(m)
+ const double Pc_max,
+ bool has_regularized)
+ : CapillaryPressureSaturation(Sr, Sg_r, Smax, Pc_max),
+ _pb(pb),
+ _m(m),
+ _has_regularized(has_regularized)
{
}
@@ -76,8 +92,24 @@ public:
double getdPcdS(const double saturation) const override;
private:
- const double _pb; ///< Entry pressure.
- const double _m; ///< Exponent m, m in [0,1]. n=1/(1-m).
+ const double _pb; ///< Entry pressure.
+ const double _m; ///< Exponent m, m in [0,1]. n=1/(1-m).
+ const bool _has_regularized; /// using regularized van Genuchten model
+ const double _xi = 1e-5; /// parameter in regularized van Genuchten model
+
+private:
+ /// Regularized van Genuchten capillary pressure-saturation Model
+ double getPcBarvGSg(double Sg) const;
+ /// Regularized van Genuchten capillary pressure-saturation Model
+ double getSBar(double Sg) const;
+ /// van Genuchten capillary pressure-saturation Model
+ double getPcvGSg(double Sg) const;
+ /// derivative dPCdS based on regularized van Genuchten capillary
+ /// pressure-saturation Model
+ double getdPcdSvGBar(double Sg) const;
+ /// derivative dPCdS based on standard van Genuchten capillary
+ /// pressure-saturation Model
+ double getdPcdSvG(const double Sg) const;
};
} // end namespace
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index f2464fa7d05b89f56ebb8ef5f20f0c3f77cb3ad0..92b659f49364e3149060bc529e5c0faa683673b7 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -21,6 +21,7 @@ APPEND_SOURCE_FILES(SOURCES HydroMechanics)
add_subdirectory(TwoPhaseFlowWithPP)
APPEND_SOURCE_FILES(SOURCES TwoPhaseFlowWithPP)
+APPEND_SOURCE_FILES(SOURCES TwoPhaseFlowWithPrho)
add_subdirectory(SmallDeformation)
APPEND_SOURCE_FILES(SOURCES SmallDeformation)
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowPrhoMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowPrhoMaterialProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cd5174a7e05b40d51367e997601dd6dac0d48a21
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowPrhoMaterialProperties.cpp
@@ -0,0 +1,138 @@
+/**
+ * \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 "CreateTwoPhaseFlowPrhoMaterialProperties.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 "TwoPhaseFlowWithPrhoMaterialProperties.h"
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPrho
+{
+std::unique_ptr<TwoPhaseFlowWithPrhoMaterialProperties>
+createTwoPhaseFlowPrhoMaterialProperties(
+ 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_PRHO__material_property__fluid}
+ auto const& fluid_config = config.getConfigSubtree("fluid");
+
+ // Get fluid properties
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__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_PRHO__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_PRHO__material_property__liquid_viscosity}
+ auto const& mu_conf = fluid_config.getConfigSubtree("liquid_viscosity");
+ auto _viscosity = MaterialLib::Fluid::createViscosityModel(mu_conf);
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__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_PRHO__material_property__porous_medium}
+ auto const& poro_config = config.getConfigSubtree("porous_medium");
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__material_property__porous_medium__porous_medium}
+ for (auto const& conf : poro_config.getConfigSubtreeList("porous_medium"))
+ {
+ //! \ogs_file_attr{prj__processes__process__TWOPHASE_FLOW_PRHO__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_PRHO__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_PRHO__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_PRHO__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_PRHO__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_PRHO__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_PRHO__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_PRHO__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<TwoPhaseFlowWithPrhoMaterialProperties>{
+ new TwoPhaseFlowWithPrhoMaterialProperties{
+ 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), std::move(_capillary_pressure_models),
+ std::move(_relative_permeability_models)}};
+}
+
+} // end namespace
+} // end namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowPrhoMaterialProperties.h b/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowPrhoMaterialProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..738a99ce5a40cb9f2fde8db82428d4299c24903f
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowPrhoMaterialProperties.h
@@ -0,0 +1,29 @@
+/**
+ * \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 "ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h"
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPrho
+{
+std::unique_ptr<TwoPhaseFlowWithPrhoMaterialProperties>
+createTwoPhaseFlowPrhoMaterialProperties(
+ BaseLib::ConfigTree const& config,
+ boost::optional<MeshLib::PropertyVector<int> const&> material_ids);
+
+} // end namespace
+} // end namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowWithPrhoProcess.cpp b/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowWithPrhoProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a4267b7dd7aa53de75b5704fe745149401faeb90
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowWithPrhoProcess.cpp
@@ -0,0 +1,119 @@
+/**
+ * \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 "CreateTwoPhaseFlowWithPrhoProcess.h"
+#include <cassert>
+#include "BaseLib/Functional.h"
+#include "MeshLib/MeshGenerators/MeshGenerator.h"
+#include "ProcessLib/Parameter/ConstantParameter.h"
+#include "ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowPrhoMaterialProperties.h"
+#include "ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h"
+#include "ProcessLib/Utils/ParseSecondaryVariables.h"
+#include "ProcessLib/Utils/ProcessUtils.h"
+#include "TwoPhaseFlowWithPrhoProcess.h"
+#include "TwoPhaseFlowWithPrhoProcessData.h"
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPrho
+{
+std::unique_ptr<Process> createTwoPhaseFlowWithPrhoProcess(
+ 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{prj__processes__process__type}
+ config.checkConfigParameter("type", "TWOPHASE_FLOW_PRHO");
+
+ DBUG("Create TwoPhaseFlowProcess with Prho model.");
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__process_variables}
+ auto const pv_config = config.getConfigSubtree("process_variables");
+
+ auto process_variables = findProcessVariables(
+ variables, pv_config,
+ {//! \ogs_file_param_special{prj__processes__process__TWOPHASE_FLOW_PRHO__process_variables__liquid_pressure}
+ "liquid_pressure",
+ //! \ogs_file_param_special{prj__processes__process__TWOPHASE_FLOW_PRHO__process_variables__overall_mass_density}
+ "overall_mass_density"});
+
+ SecondaryVariableCollection secondary_variables;
+
+ NumLib::NamedFunctionCaller named_function_caller(
+ {"TwoPhaseFlow_pressure"});
+
+ ProcessLib::parseSecondaryVariables(config, secondary_variables,
+ named_function_caller);
+ // Specific body force
+ std::vector<double> const b =
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__specific_body_force}
+ config.getConfigParameter<std::vector<double>>("specific_body_force");
+ assert(b.size() > 0 && b.size() < 4);
+ 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_PRHO__mass_lumping}
+ auto const mass_lumping = config.getConfigParameter<bool>("mass_lumping");
+ // diffusion coeff
+ auto& diff_coeff_b = findParameter<double>(
+ config,
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__diffusion_coeff_component_b}
+ "diffusion_coeff_component_b", parameters, 1);
+ auto& diff_coeff_a = findParameter<double>(
+ config,
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__diffusion_coeff_component_a}
+ "diffusion_coeff_component_a", parameters, 1);
+ auto& temperature = findParameter<double>(
+ config,
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__temperature}
+ "temperature", parameters, 1);
+
+ //! \ogs_file_param{prj__processes__process__TWOPHASE_FLOW_PRHO__material_property}
+ auto const& mat_config = config.getConfigSubtree("material_property");
+
+ auto const& mat_ids =
+ mesh.getProperties().getPropertyVector<int>("MaterialIDs");
+
+ std::unique_ptr<ProcessLib::TwoPhaseFlowWithPrho::
+ TwoPhaseFlowWithPrhoMaterialProperties>
+ material = nullptr;
+
+ boost::optional<MeshLib::PropertyVector<int> const&> material_ids;
+ if (mat_ids != nullptr)
+ {
+ INFO("The twophase flow is in heterogeneous porous media.");
+ material_ids = *mat_ids;
+ }
+ else
+ {
+ INFO("The twophase flow is in homogeneous porous media.");
+ }
+
+ material = ProcessLib::TwoPhaseFlowWithPrho::
+ createTwoPhaseFlowPrhoMaterialProperties(mat_config, material_ids);
+
+ TwoPhaseFlowWithPrhoProcessData process_data{
+ specific_body_force, has_gravity, mass_lumping, diff_coeff_b,
+ diff_coeff_a, temperature, std::move(material)};
+
+ return std::unique_ptr<Process>{new TwoPhaseFlowWithPrhoProcess{
+ 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/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowWithPrhoProcess.h b/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowWithPrhoProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..5bd873098de3fa3a0ff5faf8d973638d4a821ba2
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/CreateTwoPhaseFlowWithPrhoProcess.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 <memory>
+#include "ProcessLib/Process.h"
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPrho
+{
+std::unique_ptr<Process> createTwoPhaseFlowWithPrhoProcess(
+ 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
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..7211da2309c26958eb80beae0690542040541750
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler-impl.h
@@ -0,0 +1,232 @@
+/**
+ * \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 "TwoPhaseFlowWithPrhoLocalAssembler.h"
+
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "NumLib/Function/Interpolation.h"
+#include "TwoPhaseFlowWithPrhoProcessData.h"
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPrho
+{
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+void TwoPhaseFlowWithPrhoLocalAssembler<
+ 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);
+
+ auto Mgp =
+ local_M.template block<nonwet_pressure_size, nonwet_pressure_size>(
+ nonwet_pressure_matrix_index, nonwet_pressure_matrix_index);
+ auto Mgx = 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 Mlx = 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 Kgx = 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 Klx = 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());
+ 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));
+
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ auto const& sm = _shape_matrices[ip];
+
+ double pl_int_pt = 0.;
+ double totalrho_int_pt =
+ 0.; // total mass density of the light component
+ NumLib::shapeFunctionInterpolate(local_x, sm.N, pl_int_pt,
+ totalrho_int_pt);
+
+ const double temperature = _process_data._temperature(t, pos)[0];
+
+ double const rho_gas =
+ _process_data._material->getGasDensity(pl_int_pt, temperature);
+ double const rho_h2o =
+ _process_data._material->getLiquidDensity(pl_int_pt, temperature);
+
+ double& Sw = _ip_data[ip].sw;
+ /// Here only consider one component in gas phase
+ double const X_h2_nonwet = 1.0;
+ double& rho_h2_wet = _ip_data[ip].rho_m;
+ double& dSwdP = _ip_data[ip].dsw_dpg;
+ double& dSwdrho = _ip_data[ip].dsw_drho;
+ double& drhoh2wet = _ip_data[ip].drhom_dpg;
+ double& drhoh2wet_drho = _ip_data[ip].drhom_drho;
+ if (!_ip_data[ip].mat_property.computeConstitutiveRelation(
+ t,
+ pos,
+ material_id,
+ pl_int_pt,
+ totalrho_int_pt,
+ temperature,
+ Sw,
+ rho_h2_wet,
+ dSwdP,
+ dSwdrho,
+ drhoh2wet,
+ drhoh2wet_drho))
+ OGS_FATAL("Computation of local constitutive relation failed.");
+ double const pc = _process_data._material->getCapillaryPressure(
+ material_id, t, pos, pl_int_pt, temperature, Sw);
+
+ double const rho_wet = rho_h2o + rho_h2_wet;
+ _saturation[ip] = Sw;
+ _pressure_nonwetting[ip] = pl_int_pt + pc;
+
+ // Assemble M matrix
+ // nonwetting
+ double dPC_dSw =
+ _process_data._material->getCapillaryPressureDerivative(
+ material_id, t, pos, pl_int_pt, temperature, Sw);
+
+ double const porosity = _process_data._material->getPorosity(
+ material_id, t, pos, pl_int_pt, temperature, 0);
+
+ Mgx.noalias() += porosity * _ip_data[ip].massOperator;
+
+ Mlp.noalias() += porosity * rho_h2o * dSwdP * _ip_data[ip].massOperator;
+
+ Mlx.noalias() +=
+ porosity * (1 + dSwdrho * rho_h2o) * _ip_data[ip].massOperator;
+ double const k_rel_gas =
+ _process_data._material->getNonwetRelativePermeability(
+ t, pos, _pressure_nonwetting[ip], temperature, Sw);
+ double const mu_gas = _process_data._material->getGasViscosity(
+ _pressure_nonwetting[ip], temperature);
+ double const lambda_gas = k_rel_gas / mu_gas;
+ double const diffusion_coeff_component_h2 =
+ _process_data._diffusion_coeff_component_b(t, pos)[0];
+
+ // wet
+ double const k_rel_wet =
+ _process_data._material->getWetRelativePermeability(
+ t, pos, pl_int_pt, temperature, Sw);
+ double const mu_wet =
+ _process_data._material->getLiquidViscosity(pl_int_pt, temperature);
+ double const lambda_wet = k_rel_wet / mu_wet;
+
+ laplace_operator.noalias() = sm.dNdx.transpose() * permeability *
+ sm.dNdx * _ip_data[ip].integration_weight;
+
+ Kgp.noalias() +=
+ (rho_gas * X_h2_nonwet * lambda_gas * (1 + dPC_dSw * dSwdP) +
+ rho_h2_wet * lambda_wet) *
+ laplace_operator +
+ (Sw * porosity * diffusion_coeff_component_h2 *
+ (rho_h2o / rho_wet) * drhoh2wet) *
+ _ip_data[ip].diffusionOperator;
+ Kgx.noalias() +=
+ (rho_gas * X_h2_nonwet * lambda_gas * dPC_dSw * dSwdrho) *
+ laplace_operator +
+ (Sw * porosity * diffusion_coeff_component_h2 *
+ (rho_h2o / rho_wet) * drhoh2wet_drho) *
+ _ip_data[ip].diffusionOperator;
+ Klp.noalias() += (rho_gas * lambda_gas * (1 + dPC_dSw * dSwdP) +
+ rho_wet * lambda_wet) *
+ laplace_operator;
+
+ Klx.noalias() +=
+ (rho_gas * lambda_gas * dPC_dSw * dSwdrho) * laplace_operator;
+
+ if (_process_data._has_gravity)
+ {
+ auto const& b = _process_data._specific_body_force;
+ Bg.noalias() += (rho_gas * rho_gas * lambda_gas +
+ rho_h2_wet * rho_wet * lambda_wet) *
+ sm.dNdx.transpose() * permeability * b *
+ _ip_data[ip].integration_weight;
+ Bl.noalias() += (rho_wet * lambda_wet * rho_wet +
+ rho_gas * rho_gas * lambda_gas) *
+ sm.dNdx.transpose() * permeability * b *
+ _ip_data[ip].integration_weight;
+
+ } // end of has gravity
+ }
+ if (_process_data._has_mass_lumping)
+ {
+ for (unsigned row = 0; row < Mgp.cols(); row++)
+ {
+ for (unsigned column = 0; column < Mgp.cols(); column++)
+ {
+ if (row != column)
+ {
+ Mgx(row, row) += Mgx(row, column);
+ Mgx(row, column) = 0.0;
+ Mgp(row, row) += Mgp(row, column);
+ Mgp(row, column) = 0.0;
+ Mlx(row, row) += Mlx(row, column);
+ Mlx(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
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h
new file mode 100644
index 0000000000000000000000000000000000000000..884476fd8079afc042e1d8d532a3aec643611d7b
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h
@@ -0,0 +1,179 @@
+/**
+ * \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 <vector>
+#include "MaterialLib/PhysicalConstant.h"
+#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 "TwoPhaseFlowWithPrhoProcessData.h"
+template <typename NodalMatrixType>
+struct IntegrationPointData final
+{
+ explicit IntegrationPointData(
+ ProcessLib::TwoPhaseFlowWithPrho::
+ TwoPhaseFlowWithPrhoMaterialProperties& material_property_)
+ : mat_property(material_property_),
+ sw(1.0),
+ rho_m(0.0),
+ dsw_dpg(0.0),
+ dsw_drho(0.0),
+ drhom_dpg(0.0),
+ drhom_drho(0.0)
+ {
+ }
+ ProcessLib::TwoPhaseFlowWithPrho::TwoPhaseFlowWithPrhoMaterialProperties&
+ mat_property;
+ double sw;
+ double rho_m;
+ double dsw_dpg;
+ double dsw_drho;
+ double drhom_dpg;
+ double drhom_drho;
+ double pressure_nonwetting;
+
+ double integration_weight;
+ NodalMatrixType massOperator;
+ NodalMatrixType diffusionOperator;
+};
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPrho
+{
+const unsigned NUM_NODAL_DOF = 2;
+
+class TwoPhaseFlowWithPrhoLocalAssemblerInterface
+ : public ProcessLib::LocalAssemblerInterface,
+ public NumLib::ExtrapolatableElement
+{
+public:
+ virtual std::vector<double> const& getIntPtSaturation(
+ std::vector<double>& /*cache*/) const = 0;
+
+ virtual std::vector<double> const& getIntPtNonWettingPressure(
+ std::vector<double>& /*cache*/) const = 0;
+};
+
+template <typename ShapeFunction, typename IntegrationMethod,
+ unsigned GlobalDim>
+class TwoPhaseFlowWithPrhoLocalAssembler
+ : public TwoPhaseFlowWithPrhoLocalAssemblerInterface
+{
+ 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:
+ TwoPhaseFlowWithPrhoLocalAssembler(
+ MeshLib::Element const& element,
+ std::size_t const /*local_matrix_size*/,
+ bool const is_axially_symmetric,
+ unsigned const integration_order,
+ TwoPhaseFlowWithPrhoProcessData 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_nonwetting(
+ std::vector<double>(_integration_method.getNumberOfPoints()))
+ {
+ unsigned const n_integration_points =
+ _integration_method.getNumberOfPoints();
+ _ip_data.reserve(n_integration_points);
+ for (unsigned ip = 0; ip < n_integration_points; ip++)
+ {
+ _ip_data.emplace_back(*_process_data._material);
+ auto const& sm = _shape_matrices[ip];
+ _ip_data[ip].integration_weight =
+ sm.integralMeasure * sm.detJ *
+ _integration_method.getWeightedPoint(ip).getWeight();
+ _ip_data[ip].massOperator.setZero(ShapeFunction::NPOINTS,
+ ShapeFunction::NPOINTS);
+ _ip_data[ip].diffusionOperator.setZero(ShapeFunction::NPOINTS,
+ ShapeFunction::NPOINTS);
+ _ip_data[ip].massOperator.noalias() =
+ sm.N.transpose() * sm.N * _ip_data[ip].integration_weight;
+ _ip_data[ip].diffusionOperator.noalias() =
+ sm.dNdx.transpose() * sm.dNdx * _ip_data[ip].integration_weight;
+ }
+ }
+
+ 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& getIntPtNonWettingPressure(
+ std::vector<double>& /*cache*/) const override
+ {
+ assert(_pressure_nonwetting.size() > 0);
+ return _pressure_nonwetting;
+ }
+
+private:
+ MeshLib::Element const& _element;
+
+ IntegrationMethod const _integration_method;
+ std::vector<ShapeMatrices> _shape_matrices;
+
+ TwoPhaseFlowWithPrhoProcessData const& _process_data;
+ std::vector<IntegrationPointData<NodalMatrixType>> _ip_data;
+
+ std::vector<double> _saturation; /// used for secondary variable output
+ std::vector<double> _pressure_nonwetting;
+ 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 "TwoPhaseFlowWithPrhoLocalAssembler-impl.h"
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a719eff4603f6d45b5f16c98b526331a411e0ac8
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
@@ -0,0 +1,384 @@
+/**
+* \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 "TwoPhaseFlowWithPrhoMaterialProperties.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 "NumLib/NewtonRaphson.h"
+#include "ProcessLib/Parameter/Parameter.h"
+#include "ProcessLib/Parameter/SpatialPosition.h"
+
+using MaterialLib::PhysicalConstant::MolarMass::H2;
+using MaterialLib::PhysicalConstant::IdealGasConstant;
+using MaterialLib::PhysicalConstant::HenryConstant::HenryConstantH2;
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPrho
+{
+TwoPhaseFlowWithPrhoMaterialProperties::TwoPhaseFlowWithPrhoMaterialProperties(
+ boost::optional<MeshLib::PropertyVector<int> const&> 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,
+ 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)),
+ _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)),
+ _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 P-RHO model.");
+}
+
+int TwoPhaseFlowWithPrhoMaterialProperties::getMaterialID(
+ const std::size_t element_id)
+{
+ if (!_material_ids)
+ {
+ return 0;
+ }
+
+ assert(element_id < _material_ids->size());
+ return (*_material_ids)[element_id];
+}
+
+double TwoPhaseFlowWithPrhoMaterialProperties::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 TwoPhaseFlowWithPrhoMaterialProperties::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 TwoPhaseFlowWithPrhoMaterialProperties::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::p)] = p;
+
+ return _gas_density->getdValue(vars,
+ MaterialLib::Fluid::PropertyVariableType::p);
+}
+double TwoPhaseFlowWithPrhoMaterialProperties::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 _viscosity->getValue(vars);
+}
+
+double TwoPhaseFlowWithPrhoMaterialProperties::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& TwoPhaseFlowWithPrhoMaterialProperties::getPermeability(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const int /*dim*/) const
+{
+ return _intrinsic_permeability_models[material_id];
+}
+
+double TwoPhaseFlowWithPrhoMaterialProperties::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 TwoPhaseFlowWithPrhoMaterialProperties::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 TwoPhaseFlowWithPrhoMaterialProperties::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 TwoPhaseFlowWithPrhoMaterialProperties::getCapillaryPressure(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+ const double /*T*/, const double saturation) const
+{
+ return _capillary_pressure_models[material_id]->getCapillaryPressure(
+ saturation);
+}
+
+double TwoPhaseFlowWithPrhoMaterialProperties::getCapillaryPressureDerivative(
+ const int material_id, const double /*t*/,
+ const ProcessLib::SpatialPosition& /*pos*/, const double /*p*/,
+ const double /*T*/, const double saturation) const
+{
+ return _capillary_pressure_models[material_id]->getdPcdS(saturation);
+}
+
+bool TwoPhaseFlowWithPrhoMaterialProperties::computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ int const material_id,
+ double const pg,
+ double const X,
+ double const T,
+ double& Sw,
+ double& X_m,
+ double& dsw_dpg,
+ double& dsw_dX,
+ double& dxm_dpg,
+ double& dxm_dX)
+{
+ { // Local Newton solver
+ using LocalJacobianMatrix =
+ Eigen::Matrix<double, 2, 2, Eigen::RowMajor>;
+ using LocalResidualVector = Eigen::Matrix<double, 2, 1>;
+ using LocalUnknownVector = Eigen::Matrix<double, 2, 1>;
+ LocalJacobianMatrix J_loc;
+
+ Eigen::PartialPivLU<LocalJacobianMatrix> linear_solver(2);
+ auto const update_residual = [&](LocalResidualVector& residual) {
+ calculateResidual(material_id, pg, X, T, Sw, X_m, residual);
+ };
+
+ auto const update_jacobian = [&](LocalJacobianMatrix& jacobian) {
+ calculateJacobian(material_id, t, x, pg, X, T, jacobian, Sw,
+ X_m); // for solution dependent Jacobians
+ };
+
+ auto const update_solution = [&](LocalUnknownVector const& increment) {
+ // increment solution vectors
+ Sw += increment[0];
+ X_m += increment[1];
+ };
+
+ // TODO Make the following choice of maximum iterations and convergence
+ // criteria available from the input file configuration:
+ const int maximum_iterations(20);
+ const double tolerance(1.e-14);
+
+ auto newton_solver = NumLib::NewtonRaphson<
+ decltype(linear_solver), LocalJacobianMatrix,
+ decltype(update_jacobian), LocalResidualVector,
+ decltype(update_residual), decltype(update_solution)>(
+ linear_solver, update_jacobian, update_residual, update_solution,
+ maximum_iterations, tolerance);
+
+ auto const success_iterations = newton_solver.solve(J_loc);
+
+ if (!success_iterations)
+ return false;
+ }
+ dsw_dpg = calculatedSwdP(pg, Sw, X_m, T, material_id);
+ dsw_dX = calculatedSwdX(pg, X, Sw, X_m, T, material_id);
+ dxm_dpg = calculatedXmdP(pg, Sw, X_m, dsw_dpg, material_id);
+ dxm_dX = calculatedXmdX(pg, Sw, X_m, dsw_dX, material_id);
+ return true;
+}
+void TwoPhaseFlowWithPrhoMaterialProperties::calculateResidual(
+ const int material_id, double const pl, double const X, double const T,
+ double Sw, double rho_h2_wet, ResidualVector& res)
+{
+ const double pg =
+ pl + _capillary_pressure_models[material_id]->getCapillaryPressure(Sw);
+ const double rho_h2_nonwet = pg * H2 / IdealGasConstant / T;
+
+ // calculating residual
+ res(0) = calculateEquilibiumRhoWetLight(pg, Sw, rho_h2_wet);
+ res(1) = calculateSaturation(pl, X, Sw, rho_h2_wet, rho_h2_nonwet, T);
+}
+
+void TwoPhaseFlowWithPrhoMaterialProperties::calculateJacobian(
+ const int material_id, double const /*t*/,
+ ProcessLib::SpatialPosition const& /*x*/, double const pl,
+ double const /*X*/, double const T, JacobianMatrix& Jac, double Sw,
+ double rho_h2_wet)
+{
+ const double pg =
+ pl + _capillary_pressure_models[material_id]->getCapillaryPressure(Sw);
+ const double rho_h2_nonwet = pg * H2 / IdealGasConstant / T;
+ double const rho_equili_h2_wet = pg * HenryConstantH2 * H2;
+ double const dPC_dSw =
+ _capillary_pressure_models[material_id]->getdPcdS(Sw);
+ double const drhoh2wet_dpg = HenryConstantH2 * H2;
+ Jac.setZero();
+ if ((1 - Sw) < (rho_equili_h2_wet - rho_h2_wet))
+ {
+ Jac(0, 0) = -1;
+ }
+ else
+ {
+ Jac(0, 0) = drhoh2wet_dpg * dPC_dSw;
+ Jac(0, 1) = -1;
+ }
+
+ Jac(1, 0) = rho_h2_nonwet - rho_h2_wet;
+ Jac(1, 1) = -Sw;
+}
+
+/** Complementary condition 1
+* for calculating molar fraction of light component in the liquid phase
+*/
+double TwoPhaseFlowWithPrhoMaterialProperties::calculateEquilibiumRhoWetLight(
+ double const pg, double const Sw, double const rho_wet_h2) const
+{
+ double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
+ return std::min(1 - Sw, rho_equilibrium_wet_h2 - rho_wet_h2);
+}
+
+/** Complementary condition 2
+* for calculating the saturation
+*/
+double TwoPhaseFlowWithPrhoMaterialProperties::calculateSaturation(
+ double /*PL*/, double X, double Sw, double rho_wet_h2, double rho_nonwet_h2,
+ double /*T*/) const
+{
+ return X - (Sw * rho_wet_h2 + (1 - Sw) * rho_nonwet_h2);
+}
+
+/**
+* Calculate the derivatives using the analytical way
+*/
+double TwoPhaseFlowWithPrhoMaterialProperties::calculatedSwdP(
+ double pl, double S, double rho_wet_h2, double const T,
+ int current_material_id) const
+{
+ const double pg =
+ pl +
+ _capillary_pressure_models[current_material_id]->getCapillaryPressure(
+ S);
+ double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
+ if ((1 - S) < (rho_equilibrium_wet_h2 - rho_wet_h2))
+ {
+ return 0.0;
+ }
+ double const drhoh2wet_dpg = HenryConstantH2 * H2;
+ double const drhoh2nonwet_dpg = H2 / IdealGasConstant / T;
+ double const alpha =
+ ((drhoh2nonwet_dpg - drhoh2wet_dpg) * (1 - S) + drhoh2wet_dpg);
+ double const beta = (drhoh2nonwet_dpg - drhoh2wet_dpg) *
+ pg; // NOTE here should be PG^h, but we ignore vapor
+ double const dPC_dSw =
+ _capillary_pressure_models[current_material_id]->getdPcdS(S);
+ return alpha / (beta - alpha * dPC_dSw);
+}
+/**
+* Calculate the derivatives using the analytical way
+*/
+double TwoPhaseFlowWithPrhoMaterialProperties::calculatedSwdX(
+ double const pl, const double /*X*/, const double S,
+ const double rho_wet_h2, double const T, int current_material_id) const
+{
+ const double pg =
+ pl +
+ _capillary_pressure_models[current_material_id]->getCapillaryPressure(
+ S);
+ double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
+ if ((1 - S) < (rho_equilibrium_wet_h2 - rho_wet_h2))
+ {
+ return 0.0;
+ }
+ double const drhoh2wet_dpg = HenryConstantH2 * H2;
+ double const drhoh2nonwet_dpg = H2 / IdealGasConstant / T;
+ double const alpha =
+ ((drhoh2nonwet_dpg - drhoh2wet_dpg) * (1 - S) + drhoh2wet_dpg);
+ double const beta = (drhoh2nonwet_dpg - drhoh2wet_dpg) *
+ pg; // NOTE here should be PG^h, but we ignore vapor
+ double const dPC_dSw =
+ _capillary_pressure_models[current_material_id]->getdPcdS(S);
+ return -1 / (beta - alpha * dPC_dSw);
+}
+/**
+* Calculate the derivatives using the analytical way
+*/
+double TwoPhaseFlowWithPrhoMaterialProperties::calculatedXmdX(
+ double pl, double Sw, double rho_wet_h2, double dSwdX,
+ int current_material_id) const
+{
+ const double pg =
+ pl +
+ _capillary_pressure_models[current_material_id]->getCapillaryPressure(
+ Sw);
+ double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
+ double const dPC_dSw =
+ _capillary_pressure_models[current_material_id]->getdPcdS(Sw);
+ if ((1 - Sw) < (rho_equilibrium_wet_h2 - rho_wet_h2))
+ return 1.0;
+ return HenryConstantH2 * H2 * dPC_dSw * dSwdX;
+}
+/**
+* Calculate the derivatives using the analytical way
+*/
+double TwoPhaseFlowWithPrhoMaterialProperties::calculatedXmdP(
+ double pl, double Sw, double rho_wet_h2, double dSwdP,
+ int current_material_id) const
+{
+ const double pg =
+ pl +
+ _capillary_pressure_models[current_material_id]->getCapillaryPressure(
+ Sw);
+ double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
+ double const dPC_dSw =
+ _capillary_pressure_models[current_material_id]->getdPcdS(Sw);
+ if ((1 - Sw) < (rho_equilibrium_wet_h2 - rho_wet_h2))
+ return 0.0;
+ return HenryConstantH2 * H2 * (1 + dPC_dSw * dSwdP);
+}
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h
new file mode 100644
index 0000000000000000000000000000000000000000..fd232b493a986ed0c04c368764530e840dc4a0b7
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h
@@ -0,0 +1,192 @@
+/**
+* \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 <iostream>
+#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"
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+
+namespace MeshLib
+{
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
+namespace ProcessLib
+{
+class SpatialPosition;
+namespace TwoPhaseFlowWithPrho
+{
+class TwoPhaseFlowWithPrhoMaterialProperties
+{
+public:
+ using ArrayType = MaterialLib::Fluid::FluidProperty::ArrayType;
+
+ TwoPhaseFlowWithPrhoMaterialProperties(
+ boost::optional<MeshLib::PropertyVector<int> const&> 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,
+ 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 getCapillaryPressure(const int material_id, const double t,
+ const ProcessLib::SpatialPosition& pos,
+ const double p, const double T,
+ const double saturation) const;
+ double getCapillaryPressureDerivative(
+ 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 getDerivGasDensity(double const p, double const T) const;
+ bool computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x_position,
+ const int material_id,
+ double const pg,
+ double const X,
+ double const T,
+ double& Sw,
+ double& X_m,
+ double& dsw_dpg,
+ double& dsw_dX,
+ double& dxm_dpg,
+ double& dxm_dX);
+
+protected:
+ 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.
+ */
+ 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;
+
+private:
+ static int const jacobian_residual_size = 2;
+ using ResidualVector = Eigen::Matrix<double, jacobian_residual_size, 1>;
+ using JacobianMatrix =
+ Eigen::Matrix<double, jacobian_residual_size, jacobian_residual_size,
+ Eigen::RowMajor>;
+ using UnknownVector = Eigen::Matrix<double, jacobian_residual_size, 1>;
+
+private:
+ /**
+ * Calculates the residual vector.
+ */
+ void calculateResidual(const int material_id, double const pl,
+ double const X, double const T, double Sw,
+ double X_m, ResidualVector& res);
+ /**
+ * Calculates the Jacobian.
+ */
+ void calculateJacobian(const int material_id, double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const pl, double const X, double const T,
+ JacobianMatrix& Jac, double Sw, double X_m);
+ /** Complementary condition 1
+ * for calculating molar fraction of light component in the liquid phase
+ */
+ double calculateEquilibiumRhoWetLight(double const pg, double const Sw,
+ double const rho_wet_h2) const;
+ /** Complementary condition 2
+ * for calculating the saturation
+ */
+ double calculateSaturation(double /*PL*/, double X, double Sw,
+ double rho_wet_h2, double rho_nonwet_h2,
+ double /*T*/) const;
+ /**
+ * Calculate the derivatives using the analytical way
+ */
+ double calculatedSwdP(double pl, double S, double rho_wet_h2,
+ double const T, int current_material_id) const;
+ /**
+ * Calculate the derivatives using the analytical way
+ */
+ double calculatedSwdX(double const pl, const double /*X*/, const double S,
+ const double rho_wet_h2, double const T,
+ int current_material_id) const;
+ /**
+ * Calculate the derivatives using the analytical way
+ */
+ double calculatedXmdX(double pl, double Sw, double rho_wet_h2, double dSwdX,
+ int current_material_id) const;
+ /**
+ * Calculate the derivatives using the analytical way
+ */
+ double calculatedXmdP(double pl, double Sw, double rho_wet_h2, double dSwdP,
+ int current_material_id) const;
+};
+
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a07022810667e806dd9e865a1636d657a791aa3f
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.cpp
@@ -0,0 +1,108 @@
+/**
+ * \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 "TwoPhaseFlowWithPrhoProcess.h"
+
+#include <cassert>
+#include "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "MeshLib/PropertyVector.h"
+
+#include "ProcessLib/Utils/CreateLocalAssemblers.h"
+#include "TwoPhaseFlowWithPrhoLocalAssembler.h"
+
+#include "MaterialLib/PorousMedium/Porosity/Porosity.h"
+#include "MaterialLib/PorousMedium/Storage/Storage.h"
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPrho
+{
+TwoPhaseFlowWithPrhoProcess::TwoPhaseFlowWithPrhoProcess(
+ 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,
+ TwoPhaseFlowWithPrhoProcessData&& 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 Prho model.");
+}
+
+void TwoPhaseFlowWithPrhoProcess::initializeConcreteProcess(
+ NumLib::LocalToGlobalIndexMap const& dof_table,
+ MeshLib::Mesh const& mesh,
+ unsigned const integration_order)
+{
+ ProcessLib::ProcessVariable const& pv = getProcessVariables()[0];
+ ProcessLib::createLocalAssemblers<TwoPhaseFlowWithPrhoLocalAssembler>(
+ 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,
+ &TwoPhaseFlowWithPrhoLocalAssemblerInterface::getIntPtSaturation));
+
+ _secondary_variables.addSecondaryVariable(
+ "pressure_nonwetting", 1,
+ makeExtrapolator(getExtrapolator(), _local_assemblers,
+ &TwoPhaseFlowWithPrhoLocalAssemblerInterface::
+ getIntPtNonWettingPressure));
+}
+
+void TwoPhaseFlowWithPrhoProcess::assembleConcreteProcess(const double t,
+ GlobalVector const& x,
+ GlobalMatrix& M,
+ GlobalMatrix& K,
+ GlobalVector& b)
+{
+ DBUG("Assemble TwoPhaseFlowWithPrhoProcess.");
+ // 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 TwoPhaseFlowWithPrhoProcess::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 TwoPhaseFlowWithPrhoProcess.");
+
+ // 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);
+}
+void TwoPhaseFlowWithPrhoProcess::preTimestepConcreteProcess(
+ GlobalVector const& x, double const t, double const dt)
+{
+ DBUG("PreTimestep HydroMechanicsProcess.");
+
+ GlobalExecutor::executeMemberOnDereferenced(
+ &LocalAssemblerInterface::preTimestep, _local_assemblers,
+ *_local_to_global_index_map, x, t, dt);
+}
+
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
new file mode 100644
index 0000000000000000000000000000000000000000..ccf2163602fa2908caa449cd5e9213e6bc1ee8dc
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcess.h
@@ -0,0 +1,76 @@
+/**
+ * \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 "MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "ProcessLib/Process.h"
+#include "TwoPhaseFlowWithPrhoLocalAssembler.h"
+
+namespace MeshLib
+{
+class Element;
+class Mesh;
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
+namespace ProcessLib
+{
+namespace TwoPhaseFlowWithPrho
+{
+/**
+ * \brief A class to simulate the two-phase flow process with P-rho model in
+ * porous media
+ */
+class TwoPhaseFlowWithPrhoProcess final : public Process
+{
+public:
+ TwoPhaseFlowWithPrhoProcess(
+ 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,
+ TwoPhaseFlowWithPrhoProcessData&& 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;
+
+ void preTimestepConcreteProcess(GlobalVector const& x, const double t,
+ const double delta_t) override;
+
+ TwoPhaseFlowWithPrhoProcessData _process_data;
+
+ std::vector<std::unique_ptr<TwoPhaseFlowWithPrhoLocalAssemblerInterface>>
+ _local_assemblers;
+};
+
+} // end of namespace
+} // end of namespace
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcessData.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcessData.h
new file mode 100644
index 0000000000000000000000000000000000000000..460ea559dafa245ded0054318e3000485206f4b1
--- /dev/null
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoProcessData.h
@@ -0,0 +1,77 @@
+/**
+ * \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 "TwoPhaseFlowWithPrhoMaterialProperties.h"
+namespace MeshLib
+{
+class Element;
+}
+
+namespace ProcessLib
+{
+template <typename T>
+struct Parameter;
+
+namespace TwoPhaseFlowWithPrho
+{
+struct TwoPhaseFlowWithPrhoProcessData
+{
+ TwoPhaseFlowWithPrhoProcessData(
+ Eigen::VectorXd const specific_body_force_,
+ bool const has_gravity_,
+ bool const has_mass_lumping_,
+ Parameter<double> const& diffusion_coeff_component_b_,
+ Parameter<double> const& diffusion_coeff_component_a_,
+ Parameter<double> const& temperature_,
+ std::unique_ptr<TwoPhaseFlowWithPrhoMaterialProperties>&& material_)
+ : _specific_body_force(specific_body_force_),
+ _has_gravity(has_gravity_),
+ _has_mass_lumping(has_mass_lumping_),
+ _diffusion_coeff_component_b(diffusion_coeff_component_b_),
+ _diffusion_coeff_component_a(diffusion_coeff_component_a_),
+ _temperature(temperature_),
+ _material(std::move(material_))
+
+ {
+ }
+
+ TwoPhaseFlowWithPrhoProcessData(TwoPhaseFlowWithPrhoProcessData&& other)
+ : _specific_body_force(other._specific_body_force),
+ _has_gravity(other._has_gravity),
+ _has_mass_lumping(other._has_mass_lumping),
+ _diffusion_coeff_component_b(other._diffusion_coeff_component_b),
+ _diffusion_coeff_component_a(other._diffusion_coeff_component_a),
+ _temperature(other._temperature),
+ _material(std::move(other._material))
+ {
+ }
+
+ //! Copies are forbidden.
+ TwoPhaseFlowWithPrhoProcessData(TwoPhaseFlowWithPrhoProcessData const&) =
+ delete;
+
+ //! Assignments are not needed.
+ void operator=(TwoPhaseFlowWithPrhoProcessData const&) = delete;
+
+ //! Assignments are not needed.
+ void operator=(TwoPhaseFlowWithPrhoProcessData&&) = delete;
+ Eigen::VectorXd const _specific_body_force;
+
+ bool const _has_gravity;
+ bool const _has_mass_lumping;
+ Parameter<double> const& _diffusion_coeff_component_b;
+ Parameter<double> const& _diffusion_coeff_component_a;
+ Parameter<double> const& _temperature;
+ std::unique_ptr<TwoPhaseFlowWithPrhoMaterialProperties> _material;
+};
+
+} // namespace TwoPhaseFlowWithPrho
+} // namespace ProcessLib