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        2572 : PorousFlow2PhasePPTempl<is_ad>::validParams()
      19             : {
      20        2572 :   InputParameters params = PorousFlowVariableBaseTempl<is_ad>::validParams();
      21        5144 :   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        5144 :   params.addRequiredCoupledVar("phase1_porepressure",
      26             :                                "Variable that is the porepressure of phase 1 (eg, the gas phase)");
      27        5144 :   params.addRequiredParam<UserObjectName>("capillary_pressure",
      28             :                                           "Name of the UserObject defining the capillary pressure");
      29        2572 :   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        2572 :   return params;
      33           0 : }
      34             : 
      35             : template <bool is_ad>
      36        2007 : PorousFlow2PhasePPTempl<is_ad>::PorousFlow2PhasePPTempl(const InputParameters & parameters)
      37             :   : PorousFlowVariableBaseTempl<is_ad>(parameters),
      38             : 
      39        4014 :     _phase0_porepressure(_nodal_material
      40         990 :                              ? this->template coupledGenericDofValue<is_ad>("phase0_porepressure")
      41        3024 :                              : this->template coupledGenericValue<is_ad>("phase0_porepressure")),
      42        2007 :     _phase0_gradp_qp(this->template coupledGenericGradient<is_ad>("phase0_porepressure")),
      43        2007 :     _phase0_porepressure_varnum(coupled("phase0_porepressure")),
      44        2007 :     _p0var(_dictator.isPorousFlowVariable(_phase0_porepressure_varnum)
      45        2007 :                ? _dictator.porousFlowVariableNum(_phase0_porepressure_varnum)
      46             :                : 0),
      47             : 
      48        4014 :     _phase1_porepressure(_nodal_material
      49        2007 :                              ? this->template coupledGenericDofValue<is_ad>("phase1_porepressure")
      50        3024 :                              : this->template coupledGenericValue<is_ad>("phase1_porepressure")),
      51        2007 :     _phase1_gradp_qp(this->template coupledGenericGradient<is_ad>("phase1_porepressure")),
      52        2007 :     _phase1_porepressure_varnum(coupled("phase1_porepressure")),
      53        2007 :     _p1var(_dictator.isPorousFlowVariable(_phase1_porepressure_varnum)
      54        2007 :                ? _dictator.porousFlowVariableNum(_phase1_porepressure_varnum)
      55             :                : 0),
      56        4014 :     _pc_uo(this->template getUserObject<PorousFlowCapillaryPressure>("capillary_pressure"))
      57             : {
      58        2007 :   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        2007 : }
      64             : 
      65             : template <bool is_ad>
      66             : void
      67       30122 : PorousFlow2PhasePPTempl<is_ad>::initQpStatefulProperties()
      68             : {
      69       30122 :   PorousFlowVariableBaseTempl<is_ad>::initQpStatefulProperties();
      70       30122 :   buildQpPPSS();
      71       30122 : }
      72             : 
      73             : template <bool is_ad>
      74             : void
      75     2893992 : PorousFlow2PhasePPTempl<is_ad>::computeQpProperties()
      76             : {
      77             :   // size stuff correctly and prepare the derivative matrices with zeroes
      78     2893992 :   PorousFlowVariableBaseTempl<is_ad>::computeQpProperties();
      79             : 
      80     2893992 :   const auto pc = buildQpPPSS();
      81     2893992 :   const auto ds = _pc_uo.dSaturation(pc); // dS/d(pc)
      82             : 
      83     2893992 :   if (!_nodal_material)
      84             :   {
      85     1405064 :     (*_gradp_qp)[_qp][0] = _phase0_gradp_qp[_qp];
      86     1405064 :     (*_gradp_qp)[_qp][1] = _phase1_gradp_qp[_qp];
      87     1416374 :     (*_grads_qp)[_qp][0] = ds * ((*_gradp_qp)[_qp][0] - (*_gradp_qp)[_qp][1]);
      88     1416374 :     (*_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     2882682 :     if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
      96             :     {
      97     2880822 :       (*_dporepressure_dvar)[_qp][0][_p0var] = 1.0;
      98     2880822 :       if (!_nodal_material)
      99     1391894 :         (*_dgradp_qp_dgradv)[_qp][0][_p0var] = 1.0;
     100             :     }
     101             : 
     102     2882682 :     if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
     103             :     {
     104     2880822 :       (*_dporepressure_dvar)[_qp][1][_p1var] = 1.0;
     105     2880822 :       if (!_nodal_material)
     106     1391894 :         (*_dgradp_qp_dgradv)[_qp][1][_p1var] = 1.0;
     107             :     }
     108             : 
     109     2882682 :     if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
     110             :     {
     111     2880822 :       (*_dsaturation_dvar)[_qp][0][_p0var] = ds;
     112     2880822 :       (*_dsaturation_dvar)[_qp][1][_p0var] = -ds;
     113             :     }
     114             : 
     115     2882682 :     if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
     116             :     {
     117     2880822 :       (*_dsaturation_dvar)[_qp][0][_p1var] = -ds;
     118     2880822 :       (*_dsaturation_dvar)[_qp][1][_p1var] = ds;
     119             :     }
     120             : 
     121     2882682 :     if (!_nodal_material)
     122             :     {
     123     1393754 :       const auto d2s_qp = _pc_uo.d2Saturation(pc); // d^2(S_qp)/d(pc_qp)^2
     124             : 
     125     1393754 :       if (_dictator.isPorousFlowVariable(_phase0_porepressure_varnum))
     126             :       {
     127     1391894 :         (*_dgrads_qp_dgradv)[_qp][0][_p0var] = ds;
     128     1391894 :         (*_dgrads_qp_dv)[_qp][0][_p0var] = d2s_qp * _phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp];
     129     1391894 :         (*_dgrads_qp_dgradv)[_qp][1][_p0var] = -ds;
     130     1391894 :         (*_dgrads_qp_dv)[_qp][1][_p0var] = -d2s_qp * _phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp];
     131             :       }
     132             : 
     133     1393754 :       if (_dictator.isPorousFlowVariable(_phase1_porepressure_varnum))
     134             :       {
     135     1391894 :         (*_dgrads_qp_dgradv)[_qp][0][_p1var] = -ds;
     136     1391894 :         (*_dgrads_qp_dv)[_qp][0][_p1var] = -d2s_qp * _phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp];
     137     1391894 :         (*_dgrads_qp_dgradv)[_qp][1][_p1var] = ds;
     138     1391894 :         (*_dgrads_qp_dv)[_qp][1][_p1var] = d2s_qp * _phase0_gradp_qp[_qp] - _phase1_gradp_qp[_qp];
     139             :       }
     140             :     }
     141             :   }
     142     2893992 : }
     143             : 
     144             : template <bool is_ad>
     145             : GenericReal<is_ad>
     146     2924114 : PorousFlow2PhasePPTempl<is_ad>::buildQpPPSS()
     147             : {
     148     2924114 :   _porepressure[_qp][0] = _phase0_porepressure[_qp];
     149     2924114 :   _porepressure[_qp][1] = _phase1_porepressure[_qp];
     150             : 
     151     2912674 :   const auto pc = _phase0_porepressure[_qp] - _phase1_porepressure[_qp]; // this is <= 0
     152     2924114 :   const auto sat = _pc_uo.saturation(pc);
     153     2924114 :   _saturation[_qp][0] = sat;
     154     2935554 :   _saturation[_qp][1] = 1.0 - sat;
     155             : 
     156     2924114 :   return pc;
     157             : }
     158             : 
     159             : template class PorousFlow2PhasePPTempl<false>;
     160             : template class PorousFlow2PhasePPTempl<true>;