Add capillary pressure output to PorousFlowProperyAux

modules/porous_flow/doc/content/source/auxkernels/PorousFlowPropertyAux.md

@@ -1,4 +1,5 @@
 # PorousFlowPropertyAux
+
 !syntax description /AuxKernels/PorousFlowPropertyAux
 
 This `AuxKernel` provides simplified access to fluid and material properties. The
@@ -11,6 +12,7 @@ following properties are available using the `property` input parameter:
 - `viscosity`
 - `mass_fraction`
 - `relperm`
+- `capillary_pressure`
 - `enthalpy`
 - `internal_energy`
 - `secondary_concentration` (m$^{3}$(secondary species)/m$^{3}$(fluid))

modules/porous_flow/include/auxkernels/PorousFlowPropertyAux.h

@@ -81,6 +81,7 @@ class PorousFlowPropertyAux : public AuxKernel
     VISCOSITY,
     MASS_FRACTION,
     RELPERM,
+    CAPILLARY_PRESSURE,
     ENTHALPY,
     INTERNAL_ENERGY,
     SECONDARY_CONCENTRATION,
@@ -91,6 +92,12 @@ class PorousFlowPropertyAux : public AuxKernel
   /// Phase index
   const unsigned int _phase;
 
+  /// Liquid phase index
+  const unsigned int _liquid_phase;
+
+  /// Gas phase index
+  const unsigned int _gas_phase;
+
   /// Fluid component index
   const unsigned int _fluid_component;
 

modules/porous_flow/src/auxkernels/PorousFlowPropertyAux.C

@@ -19,11 +19,15 @@ validParams<PorousFlowPropertyAux>()
   params.addRequiredParam<UserObjectName>(
       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
   MooseEnum property_enum("pressure saturation temperature density viscosity mass_fraction relperm "
-                          "enthalpy internal_energy secondary_concentration mineral_concentration "
-                          "mineral_reaction_rate");
+                          "capillary_pressure enthalpy internal_energy secondary_concentration "
+                          "mineral_concentration mineral_reaction_rate");
   params.addRequiredParam<MooseEnum>(
       "property", property_enum, "The fluid property that this auxillary kernel is to calculate");
   params.addParam<unsigned int>("phase", 0, "The index of the phase this auxillary kernel acts on");
+  params.addParam<unsigned int>(
+      "liquid_phase", 0, "The index of the liquid phase (used for capillary pressure)");
+  params.addParam<unsigned int>(
+      "gas_phase", 1, "The index of the gas phase (used for capillary pressure)");
   params.addParam<unsigned int>(
       "fluid_component", 0, "The index of the fluid component this auxillary kernel acts on");
   params.addParam<unsigned int>("secondary_species", 0, "The secondary chemical species number");
@@ -39,6 +43,8 @@ PorousFlowPropertyAux::PorousFlowPropertyAux(const InputParameters & parameters)
     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
     _property_enum(getParam<MooseEnum>("property").getEnum<PropertyEnum>()),
     _phase(getParam<unsigned int>("phase")),
+    _liquid_phase(getParam<unsigned int>("liquid_phase")),
+    _gas_phase(getParam<unsigned int>("gas_phase")),
     _fluid_component(getParam<unsigned int>("fluid_component")),
     _secondary_species(getParam<unsigned int>("secondary_species")),
     _mineral_species(getParam<unsigned int>("mineral_species"))
@@ -54,6 +60,24 @@ PorousFlowPropertyAux::PorousFlowPropertyAux(const InputParameters & parameters)
                "Fluid component number entered is greater than the number of fluid components "
                "specified in the Dictator. Remember that indexing starts at 0");
 
+  // Check the parameters used to calculate capillary pressure
+  if (_property_enum == PropertyEnum::CAPILLARY_PRESSURE)
+  {
+    if (_liquid_phase >= _dictator.numPhases())
+      paramError(
+          "liquid_phase",
+          "Liquid phase number entered is greater than the number of phases specified in the "
+          "Dictator. Remember that indexing starts at 0");
+
+    if (_gas_phase >= _dictator.numPhases())
+      paramError("gas_phase",
+                 "Gas phase number entered is greater than the number of phases specified in the "
+                 "Dictator. Remember that indexing starts at 0");
+
+    if (_liquid_phase == _gas_phase)
+      paramError("liquid_phase", "Liquid phase number entered cannot be equal to gas_phase");
+  }
+
   if (_property_enum == PropertyEnum::SECONDARY_CONCENTRATION &&
       (_secondary_species >= _dictator.numAqueousEquilibrium()))
     paramError("secondary_species",
@@ -101,6 +125,10 @@ PorousFlowPropertyAux::PorousFlowPropertyAux(const InputParameters & parameters)
           &getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp");
       break;
 
+    case PropertyEnum::CAPILLARY_PRESSURE:
+      _pressure = &getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp");
+      break;
+
     case PropertyEnum::ENTHALPY:
       _enthalpy = &getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_qp");
       break;
@@ -161,6 +189,10 @@ PorousFlowPropertyAux::computeValue()
       property = (*_relative_permeability)[_qp][_phase];
       break;
 
+    case PropertyEnum::CAPILLARY_PRESSURE:
+      property = (*_pressure)[_qp][_gas_phase] - (*_pressure)[_qp][_liquid_phase];
+      break;
+
     case PropertyEnum::ENTHALPY:
       property = (*_enthalpy)[_qp][_phase];
       break;

modules/porous_flow/test/tests/aux_kernels/properties.i

@@ -38,6 +38,10 @@
     order = CONSTANT
     family = MONOMIAL
   [../]
+  [./capillary_pressure]
+    order = CONSTANT
+    family = MONOMIAL
+  [../]
   [./saturation_water]
     order = CONSTANT
     family = MONOMIAL
@@ -100,6 +104,12 @@
     phase = 1
     execute_on = timestep_end
   [../]
+  [./capillary_pressure]
+    type = PorousFlowPropertyAux
+    variable = capillary_pressure
+    property = capillary_pressure
+    execute_on = timestep_end
+  [../]
   [./saturation_water]
     type = PorousFlowPropertyAux
     variable = saturation_water