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 "PorousFlowHystereticInfo.h"
      11             : 
      12             : registerMooseObject("PorousFlowApp", PorousFlowHystereticInfo);
      13             : 
      14             : InputParameters
      15        3900 : PorousFlowHystereticInfo::validParams()
      16             : {
      17        3900 :   InputParameters params = PorousFlowHystereticCapillaryPressure::validParams();
      18        7800 :   params.addCoupledVar("pc_var",
      19             :                        0,
      20             :                        "Variable that represents capillary pressure.  Depending on info_required, "
      21             :                        "this may not be used to compute the info");
      22        7800 :   params.addRequiredCoupledVar(
      23             :       "sat_var",
      24             :       "Variable that represent liquid saturation.  This is always needed to "
      25             :       "ensure the hysteretic order is computed correctly");
      26             :   MooseEnum property_enum("pc sat sat_given_pc dS_dPc_err dPc_dS_err d2S_dPc2_err d2Pc_dS2_err",
      27        7800 :                           "pc");
      28        7800 :   params.addParam<MooseEnum>(
      29             :       "info_required",
      30             :       property_enum,
      31             :       "The type of information required.  pc: capillary pressure given the saturation (pc_var is "
      32             :       "not used in this case).  sat: given the liquid saturation, compute the capillary pressure, "
      33             :       "then invert the relationship to yield liquid saturation again (pc_var is not used in this "
      34             :       "case).  This is useful to understand the non-invertibility of the hysteretic relationships. "
      35             :       " sat_given_pc: given the capillary pressure, compute the saturation.  dS_dPc_err: relative "
      36             :       "error in d(sat)/d(pc) calculation, ie (S(Pc + fd_eps) - S(Pc - fd_eps))/(2 * eps * S'(Pc)) "
      37             :       "- 1, where S' is the coded derivative (Pc is computed from sat in this case: pc_var is not "
      38             :       "used).  This is useful for checking derivatives.  dPc_dS_err: relative error in "
      39             :       "d(pc)/d(sat) calculation, ie (Pc(S + fd_eps) - Pc(S - fd_eps)) / (2 * eps * Pc'(S)) - 1, "
      40             :       "where Pc' is the coded derative (pc_var is not used in this case).  d2S_dPc2_err: relative "
      41             :       "error in d^2(sat)/d(pc)^2 (Pc is computed from sat in this case: pc_var is not used).  "
      42             :       "d2Pc_dS2_err: relative error in d^2(pc)/d(sat)^2.");
      43        7800 :   params.addParam<Real>(
      44             :       "fd_eps",
      45        7800 :       1E-8,
      46             :       "Small quantity used in computing the finite-difference approximations to derivatives");
      47        3900 :   params.addClassDescription(
      48             :       "This Material computes capillary pressure or saturation, etc.  It is primarily "
      49             :       "of use when users desire to compute hysteretic quantities such as capillary pressure for "
      50             :       "visualisation purposes.  The "
      51             :       "result is written into PorousFlow_hysteretic_info_nodal or "
      52             :       "PorousFlow_hysteretic_info_qp (depending on the at_nodes flag).  It "
      53             :       "does not "
      54             :       "compute porepressure and should not be used in simulations that employ "
      55             :       "PorousFlow*PhaseHys* "
      56             :       "Materials.");
      57        3900 :   return params;
      58        3900 : }
      59             : 
      60        3000 : PorousFlowHystereticInfo::PorousFlowHystereticInfo(const InputParameters & parameters)
      61             :   : PorousFlowHystereticCapillaryPressure(parameters),
      62           0 :     _pc(_nodal_material ? declareProperty<Real>("PorousFlow_hysteretic_capillary_pressure_nodal")
      63        6000 :                         : declareProperty<Real>("PorousFlow_hysteretic_capillary_pressure_qp")),
      64        3000 :     _info(_nodal_material ? declareProperty<Real>("PorousFlow_hysteretic_info_nodal")
      65        6000 :                           : declareProperty<Real>("PorousFlow_hysteretic_info_qp")),
      66        6000 :     _pc_val(_nodal_material ? coupledDofValues("pc_var") : coupledValue("pc_var")),
      67        6000 :     _sat_val(_nodal_material ? coupledDofValues("sat_var") : coupledValue("sat_var")),
      68        6000 :     _fd_eps(getParam<Real>("fd_eps")),
      69        9000 :     _info_enum(getParam<MooseEnum>("info_required").getEnum<InfoTypeEnum>())
      70             : {
      71             :   // note that _num_phases must be positive, otherwise get a problem with using
      72             :   // PorousFlowHysteresisOrder, so _num_phases > 0 needn't be checked here
      73        3000 : }
      74             : 
      75             : void
      76      204750 : PorousFlowHystereticInfo::initQpStatefulProperties()
      77             : {
      78      204750 :   PorousFlowHystereticCapillaryPressure::initQpStatefulProperties();
      79      204750 :   _saturation[_qp][0] = _sat_val[_qp];
      80      204750 :   _pc[_qp] = capillaryPressureQp(_sat_val[_qp]);
      81      204750 :   computeQpInfo();
      82      204750 : }
      83             : 
      84             : void
      85      446616 : PorousFlowHystereticInfo::computeQpProperties()
      86             : {
      87      446616 :   PorousFlowHystereticCapillaryPressure::computeQpProperties();
      88      446616 :   _saturation[_qp][0] = _sat_val[_qp];
      89      446616 :   _pc[_qp] = capillaryPressureQp(_sat_val[_qp]);
      90      446616 :   computeQpInfo();
      91      446616 : }
      92             : 
      93             : void
      94      651366 : PorousFlowHystereticInfo::computeQpInfo()
      95             : {
      96      651366 :   switch (_info_enum)
      97             :   {
      98      192986 :     case InfoTypeEnum::PC:
      99      192986 :       _info[_qp] = capillaryPressureQp(_sat_val[_qp]);
     100      192986 :       break;
     101      133660 :     case InfoTypeEnum::SAT:
     102             :     {
     103      133660 :       const Real pc = capillaryPressureQp(_sat_val[_qp]);
     104      133660 :       _info[_qp] = liquidSaturationQp(pc);
     105      133660 :       break;
     106             :     }
     107           0 :     case InfoTypeEnum::SAT_GIVEN_PC:
     108           0 :       _info[_qp] = liquidSaturationQp(_pc_val[_qp]);
     109           0 :       break;
     110       81180 :     case InfoTypeEnum::DS_DPC_ERR:
     111             :     {
     112       81180 :       const Real pc = capillaryPressureQp(_sat_val[_qp]);
     113             :       const Real fd =
     114       81180 :           0.5 * (liquidSaturationQp(pc + _fd_eps) - liquidSaturationQp(pc - _fd_eps)) / _fd_eps;
     115       81180 :       _info[_qp] = relativeError(fd, dliquidSaturationQp(pc));
     116       81180 :       break;
     117             :     }
     118       81180 :     case InfoTypeEnum::DPC_DS_ERR:
     119             :     {
     120       81180 :       const Real fd = 0.5 *
     121       81180 :                       (capillaryPressureQp(_sat_val[_qp] + _fd_eps) -
     122       81180 :                        capillaryPressureQp(_sat_val[_qp] - _fd_eps)) /
     123       81180 :                       _fd_eps;
     124       81180 :       _info[_qp] = relativeError(fd, dcapillaryPressureQp(_sat_val[_qp]));
     125       81180 :       break;
     126             :     }
     127       81180 :     case InfoTypeEnum::D2S_DPC2_ERR:
     128             :     {
     129       81180 :       const Real pc = capillaryPressureQp(_sat_val[_qp]);
     130             :       const Real fd =
     131       81180 :           0.5 * (dliquidSaturationQp(pc + _fd_eps) - dliquidSaturationQp(pc - _fd_eps)) / _fd_eps;
     132       81180 :       _info[_qp] = relativeError(fd, d2liquidSaturationQp(pc));
     133       81180 :       break;
     134             :     }
     135       81180 :     case InfoTypeEnum::D2PC_DS2_ERR:
     136             :     {
     137       81180 :       const Real fd = 0.5 *
     138       81180 :                       (dcapillaryPressureQp(_sat_val[_qp] + _fd_eps) -
     139       81180 :                        dcapillaryPressureQp(_sat_val[_qp] - _fd_eps)) /
     140       81180 :                       _fd_eps;
     141       81180 :       _info[_qp] = relativeError(fd, d2capillaryPressureQp(_sat_val[_qp]));
     142       81180 :       break;
     143             :     }
     144             :   }
     145      651366 : }
     146             : 
     147             : Real
     148      324720 : PorousFlowHystereticInfo::relativeError(Real finite_difference, Real hand_coded)
     149             : {
     150      324720 :   if (finite_difference == 0.0 && hand_coded == 0.0)
     151             :     return 0.0;
     152      308074 :   else if (finite_difference > 0.0 && hand_coded == 0.0)
     153             :     return std::numeric_limits<Real>::max();
     154      308074 :   else if (finite_difference < 0.0 && hand_coded == 0.0)
     155             :     return std::numeric_limits<Real>::min();
     156      299300 :   return finite_difference / hand_coded - 1.0;
     157             : }