Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "PorousFlow2PhasePP.h"
      11             : #include "PorousFlowCapillaryPressure.h"
      12             : 
      13             : registerMooseObject("PorousFlowApp", PorousFlow2PhasePP);
      14             : registerMooseObject("PorousFlowApp", ADPorousFlow2PhasePP);
      15             : 
      16             : template <bool is_ad>
      17             : InputParameters
      18        4608 : PorousFlow2PhasePPTempl<is_ad>::validParams()
      19             : {
      20        4608 :   InputParameters params = PorousFlowVariableBaseTempl<is_ad>::validParams();
      21        9216 :   params.addRequiredCoupledVar("phase0_porepressure",
      22             :                                "Variable that is the porepressure of phase "
      23             :                                "0 (eg, the water phase).  It will be <= "
      24             :                                "phase1_porepressure.");
      25        9216 :   params.addRequiredCoupledVar("phase1_porepressure",
      26             :                                "Variable that is the porepressure of phase 1 (eg, the gas phase)");
      27        9216 :   params.addRequiredParam<UserObjectName>("capillary_pressure",
      28             :                                           "Name of the UserObject defining the capillary pressure");
      29        4608 :   params.addClassDescription("This Material calculates the 2 porepressures and the 2 saturations "
      30             :                              "in a 2-phase situation, and derivatives of these with "
      31             :                              "respect to the PorousFlowVariables");
      32        4608 :   return params;
      33           0 : }
      34             : 
      35             : template <bool is_ad>
      36        3651 : PorousFlow2PhasePPTempl<is_ad>::PorousFlow2PhasePPTempl(const InputParameters & parameters)
      37             :   : PorousFlowVariableBaseTempl<is_ad>(parameters),
      38             : 
      39        7302 :     _phase0_porepressure(_nodal_material
      40        1776 :                              ? this->template coupledGenericDofValue<is_ad>("phase0_porepressure")
      41        5526 :                              : this->template coupledGenericValue<is_ad>("phase0_porepressure")),
      42        3651 :     _phase0_gradp_qp(this->template coupledGenericGradient<is_ad>("phase0_porepressure")),
      43        3651 :     _phase0_porepressure_varnum(coupled("phase0_porepressure")),
      44        3651 :     _p0var(_dictator.isPorousFlowVariable(_phase0_porepressure_varnum)
      45        3651 :                ? _dictator.porousFlowVariableNum(_phase0_porepressure_varnum)
      46             :                : 0),
      47             : 
      48        7302 :     _phase1_porepressure(_nodal_material
      49        3651 :                              ? this->template coupledGenericDofValue<is_ad>("phase1_porepressure")
      50        5526 :                              : this->template coupledGenericValue<is_ad>("phase1_porepressure")),
      51        3651 :     _phase1_gradp_qp(this->template coupledGenericGradient<is_ad>("phase1_porepressure")),
      52        3651 :     _phase1_porepressure_varnum(coupled("phase1_porepressure")),
      53        3651 :     _p1var(_dictator.isPorousFlowVariable(_phase1_porepressure_varnum)
      54        3651 :                ? _dictator.porousFlowVariableNum(_phase1_porepressure_varnum)
      55             :                : 0),
      56        7302 :     _pc_uo(this->template getUserObject<PorousFlowCapillaryPressure>("capillary_pressure"))
      57             : {
      58        3651 :   if (_num_phases != 2)
      59           0 :     mooseError("The Dictator announces that the number of phases is ",
      60           0 :                _dictator.numPhases(),
      61             :                " whereas PorousFlow2PhasePP can only be used for 2-phase simulation.  When you "
      62             :                "have an efficient government, you have a dictatorship.");
      63        3651 : }
      64             : 
      65             : template <bool is_ad>
      66             : void
      67       47402 : PorousFlow2PhasePPTempl<is_ad>::initQpStatefulProperties()
      68             : {
      69       47402 :   PorousFlowVariableBaseTempl<is_ad>::initQpStatefulProperties();
      70       47402 :   buildQpPPSS();
      71       47402 : }
      72             : 
      73             : template <bool is_ad>
      74             : void
      75     4205133 : PorousFlow2PhasePPTempl<is_ad>::computeQpProperties()
      76             : {
      77             :   // size stuff correctly and prepare the derivative matrices with zeroes
      78     4205133 :   PorousFlowVariableBaseTempl<is_ad>::computeQpProperties();
      79             : 
      80     4205133 :   const auto pc = buildQpPPSS();
      81     4205133 :   const auto ds = _pc_uo.dSaturation(pc); // dS/d(pc)
      82             : 
      83     4205133 :   if (!_nodal_material)
      84             :   {
      85     2053471 :     (*_gradp_qp)[_qp][0] = _phase0_gradp_qp[_qp];
      86     2053471 :     (*_gradp_qp)[_qp][1] = _phase1_gradp_qp[_qp];
      87     2070361 :     (*_grads_qp)[_qp][0] = ds * ((*_gradp_qp)[_qp][0] - (*_gradp_qp)[_qp][1]);
      88     2070361 :     (*_grads_qp)[_qp][1] = -(*_grads_qp)[_qp][0];
      89             :   }
      90             : 
      91             :   // the derivatives of porepressure with respect to porepressure
      92             :   // remain fixed (at unity) throughout the simulation
      93             :   if constexpr (!is_ad)
      94             :   {
      95     4188243 :     if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
      96             :     {
      97     4185483 :       (*_dporepressure_dvar)[_qp][0][_p0var] = 1.0;
      98     4185483 :       if (!_nodal_material)
      99     2033821 :         (*_dgradp_qp_dgradv)[_qp][0][_p0var] = 1.0;
     100             :     }
     101             : 
     102     4188243 :     if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
     103             :     {
     104     4185483 :       (*_dporepressure_dvar)[_qp][1][_p1var] = 1.0;
     105     4185483 :       if (!_nodal_material)
     106     2033821 :         (*_dgradp_qp_dgradv)[_qp][1][_p1var] = 1.0;
     107             :     }
     108             : 
     109     4188243 :     if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
     110             :     {
     111     4185483 :       (*_dsaturation_dvar)[_qp][0][_p0var] = ds;
     112     4185483 :       (*_dsaturation_dvar)[_qp][1][_p0var] = -ds;
     113             :     }
     114             : 
     115     4188243 :     if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
     116             :     {
     117     4185483 :       (*_dsaturation_dvar)[_qp][0][_p1var] = -ds;
     118     4185483 :       (*_dsaturation_dvar)[_qp][1][_p1var] = ds;
     119             :     }
     120             : 
     121     4188243 :     if (!_nodal_material)
     122             :     {
     123     2036581 :       const auto d2s_qp = _pc_uo.d2Saturation(pc); // d^2(S_qp)/d(pc_qp)^2
     124             : 
     125     2036581 :       if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
     126             :       {
     127     2033821 :         (*_dgrads_qp_dgradv)[_qp][0][_p0var] = ds;
     128     2033821 :         (*_dgrads_qp_dv)[_qp][0][_p0var] = d2s_qp * _phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp];
     129     2033821 :         (*_dgrads_qp_dgradv)[_qp][1][_p0var] = -ds;
     130     2033821 :         (*_dgrads_qp_dv)[_qp][1][_p0var] = -d2s_qp * _phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp];
     131             :       }
     132             : 
     133     2036581 :       if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
     134             :       {
     135     2033821 :         (*_dgrads_qp_dgradv)[_qp][0][_p1var] = -ds;
     136     2033821 :         (*_dgrads_qp_dv)[_qp][0][_p1var] = -d2s_qp * _phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp];
     137     2033821 :         (*_dgrads_qp_dgradv)[_qp][1][_p1var] = ds;
     138     2033821 :         (*_dgrads_qp_dv)[_qp][1][_p1var] = d2s_qp * _phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp];
     139             :       }
     140             :     }
     141             :   }
     142     4205133 : }
     143             : 
     144             : template <bool is_ad>
     145             : GenericReal<is_ad>
     146     4252535 : PorousFlow2PhasePPTempl<is_ad>::buildQpPPSS()
     147             : {
     148     4252535 :   _porepressure[_qp][0] = _phase0_porepressure[_qp];
     149     4252535 :   _porepressure[_qp][1] = _phase1_porepressure[_qp];
     150             : 
     151     4252535 :   const auto pc = _phase0_porepressure[_qp] - _phase1_porepressure[_qp]; // this is <= 0
     152     4252535 :   const auto sat = _pc_uo.saturation(pc);
     153     4252535 :   _saturation[_qp][0] = sat;
     154     4269635 :   _saturation[_qp][1] = 1.0 - sat;
     155             : 
     156     4252535 :   return pc;
     157             : }
     158             : 
     159             : template class PorousFlow2PhasePPTempl<false>;
     160             : template class PorousFlow2PhasePPTempl<true>;