LCOV - code coverage report
Current view: top level - src/userobjects - PorousFlowCapillaryPressure.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 61 61 100.0 %
Date: 2017-11-21 14:47:27 Functions: 15 15 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /****************************************************************/
       2             : /* MOOSE - Multiphysics Object Oriented Simulation Environment  */
       3             : /*                                                              */
       4             : /*          All contents are licensed under LGPL V2.1           */
       5             : /*             See LICENSE for full restrictions                */
       6             : /****************************************************************/
       7             : 
       8             : #include "PorousFlowCapillaryPressure.h"
       9             : 
      10             : template <>
      11             : InputParameters
      12         193 : validParams<PorousFlowCapillaryPressure>()
      13             : {
      14         193 :   InputParameters params = validParams<GeneralUserObject>();
      15         772 :   params.addRangeCheckedParam<Real>(
      16             :       "sat_lr",
      17             :       0.0,
      18             :       "sat_lr >= 0 & sat_lr <= 1",
      19         193 :       "Liquid residual saturation.  Must be between 0 and 1. Default is 0");
      20         772 :   params.addRangeCheckedParam<Real>("pc_max",
      21             :                                     1.0e9,
      22             :                                     "pc_max > 0",
      23         193 :                                     "Maximum capillary pressure (Pa). Must be > 0. Default is 1e9");
      24         579 :   params.addParam<bool>("log_extension",
      25             :                         true,
      26             :                         "Use a logarithmic extension for low saturation to avoid capillary "
      27         193 :                         "pressure going to infinity. Default is true");
      28         386 :   params.addClassDescription("Capillary pressure base class");
      29         193 :   return params;
      30             : }
      31             : 
      32         193 : PorousFlowCapillaryPressure::PorousFlowCapillaryPressure(const InputParameters & parameters)
      33             :   : GeneralUserObject(parameters),
      34         579 :     _sat_lr(getParam<Real>("sat_lr")),
      35         193 :     _dseff_ds(1.0 / (1.0 - _sat_lr)),
      36         579 :     _log_ext(getParam<bool>("log_extension")),
      37         579 :     _pc_max(getParam<Real>("pc_max")),
      38             :     _sat_ext(0.0),
      39             :     _pc_ext(0.0),
      40             :     _slope_ext(0.0),
      41        1158 :     _log10(std::log(10.0))
      42             : {
      43         193 : }
      44             : 
      45             : void
      46         172 : PorousFlowCapillaryPressure::initialSetup()
      47             : {
      48             :   // If _log_ext = true, calculate the saturation where the the logarithmic
      49             :   // extension meets the raw capillary pressure curve.
      50         172 :   if (_log_ext)
      51             :   {
      52         136 :     _sat_ext = extensionSaturation();
      53         136 :     _pc_ext = capillaryPressure(_sat_ext);
      54         136 :     _slope_ext = (std::log10(_pc_ext) - std::log10(_pc_max)) / _sat_ext;
      55             :   }
      56         172 : }
      57             : 
      58             : Real
      59      551711 : PorousFlowCapillaryPressure::capillaryPressure(Real saturation) const
      60             : {
      61      551711 :   if (_log_ext && saturation < _sat_ext)
      62        6704 :     return capillaryPressureLogExt(saturation);
      63             :   else
      64      545007 :     return capillaryPressureCurve(saturation);
      65             : }
      66             : 
      67             : Real
      68      523656 : PorousFlowCapillaryPressure::dCapillaryPressure(Real saturation) const
      69             : {
      70      523656 :   if (_log_ext && saturation < _sat_ext)
      71        6624 :     return dCapillaryPressureLogExt(saturation);
      72             :   else
      73      517032 :     return dCapillaryPressureCurve(saturation);
      74             : }
      75             : 
      76             : Real
      77      273916 : PorousFlowCapillaryPressure::d2CapillaryPressure(Real saturation) const
      78             : {
      79      273916 :   if (_log_ext && saturation < _sat_ext)
      80        3758 :     return d2CapillaryPressureLogExt(saturation);
      81             :   else
      82      270158 :     return d2CapillaryPressureCurve(saturation);
      83             : }
      84             : 
      85             : Real
      86      209767 : PorousFlowCapillaryPressure::effectiveSaturationFromSaturation(Real saturation) const
      87             : {
      88      209767 :   return (saturation - _sat_lr) / (1.0 - _sat_lr);
      89             : }
      90             : 
      91             : Real
      92        6704 : PorousFlowCapillaryPressure::capillaryPressureLogExt(Real saturation) const
      93             : {
      94        6704 :   return _pc_ext * std::pow(10.0, _slope_ext * (saturation - _sat_ext));
      95             : }
      96             : 
      97             : Real
      98        6624 : PorousFlowCapillaryPressure::dCapillaryPressureLogExt(Real saturation) const
      99             : {
     100        6624 :   return _pc_ext * _slope_ext * _log10 * std::pow(10.0, _slope_ext * (saturation - _sat_ext));
     101             : }
     102             : 
     103             : Real
     104        3758 : PorousFlowCapillaryPressure::d2CapillaryPressureLogExt(Real saturation) const
     105             : {
     106        3758 :   return _pc_ext * _slope_ext * _slope_ext * _log10 * _log10 *
     107        3758 :          std::pow(10.0, _slope_ext * (saturation - _sat_ext));
     108             : }
     109             : 
     110             : Real
     111         136 : PorousFlowCapillaryPressure::extensionSaturation() const
     112             : {
     113             :   // Initial guess for saturation where extension matches curve
     114         136 :   Real sat = _sat_lr + 0.01;
     115             : 
     116             :   // Calculate the saturation where the extension matches the derivative of the
     117             :   // raw capillary pressure curve
     118             :   const unsigned int max_its = 20;
     119             :   const Real nr_tol = 1.0e-8;
     120             :   unsigned int iter = 0;
     121             : 
     122         728 :   while (std::abs(interceptFunction(sat)) > nr_tol)
     123             :   {
     124         228 :     sat = sat - interceptFunction(sat) / interceptFunctionDeriv(sat);
     125             : 
     126         228 :     iter++;
     127         228 :     if (iter > max_its)
     128             :       break;
     129             :   }
     130             : 
     131         136 :   return sat;
     132             : }
     133             : 
     134             : Real
     135         592 : PorousFlowCapillaryPressure::interceptFunction(Real saturation) const
     136             : {
     137         592 :   Real pc = capillaryPressureCurve(saturation);
     138         592 :   Real dpc = dCapillaryPressureCurve(saturation);
     139             : 
     140         592 :   return std::log10(pc) - saturation * dpc / (_log10 * pc) - std::log10(_pc_max);
     141             : }
     142             : 
     143             : Real
     144         228 : PorousFlowCapillaryPressure::interceptFunctionDeriv(Real saturation) const
     145             : {
     146         228 :   Real pc = capillaryPressureCurve(saturation);
     147         228 :   Real dpc = dCapillaryPressureCurve(saturation);
     148         228 :   Real d2pc = d2CapillaryPressureCurve(saturation);
     149             : 
     150         228 :   return saturation * (dpc * dpc / pc - d2pc) / (_log10 * pc);
     151        2499 : }

Generated by: LCOV version 1.11