Line data    Source code

       1              : // Copyright 2024, UChicago Argonne, LLC
       2              : // All Rights Reserved
       3              : // Software Name: NEML2 -- the New Engineering material Model Library, version 2
       4              : // By: Argonne National Laboratory
       5              : // OPEN SOURCE LICENSE (MIT)
       6              : //
       7              : // Permission is hereby granted, free of charge, to any person obtaining a copy
       8              : // of this software and associated documentation files (the "Software"), to deal
       9              : // in the Software without restriction, including without limitation the rights
      10              : // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
      11              : // copies of the Software, and to permit persons to whom the Software is
      12              : // furnished to do so, subject to the following conditions:
      13              : //
      14              : // The above copyright notice and this permission notice shall be included in
      15              : // all copies or substantial portions of the Software.
      16              : //
      17              : // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
      18              : // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      19              : // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
      20              : // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      21              : // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
      22              : // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
      23              : // THE SOFTWARE.
      24              : 
      25              : #include "neml2/models/porous_flow/CapillaryPressure.h"
      26              : #include "neml2/tensors/functions/log10.h"
      27              : #include "neml2/tensors/functions/where.h"
      28              : 
      29              : namespace neml2
      30              : {
      31              : OptionSet
      32            4 : CapillaryPressure::expected_options()
      33              : {
      34            4 :   OptionSet options = Model::expected_options();
      35            4 :   options.doc() = "Relate the porous flow capillary pressure to the effective saturation";
      36              : 
      37           12 :   options.set_input("effective_saturation") = VariableName(STATE, "S");
      38            4 :   options.set("effective_saturation").doc() = "The effective saturation";
      39              : 
      40           12 :   options.set_output("capillary_pressure") = VariableName(STATE, "Pc");
      41            8 :   options.set("capillary_pressure").doc() = "Capillary pressure.";
      42              : 
      43            8 :   options.set<bool>("log_extension") = false;
      44            8 :   options.set("log_extension").doc() = "Whether to apply logarithmic extension";
      45              : 
      46            8 :   options.set<double>("transition_saturation") = 0.0;
      47            8 :   options.set("transition_saturation").doc() = "The transistion value of the effective saturation "
      48            4 :                                                "below which to apply the logarithmic extension";
      49              : 
      50            4 :   return options;
      51            0 : }
      52              : 
      53            2 : CapillaryPressure::CapillaryPressure(const OptionSet & options)
      54              :   : Model(options),
      55            2 :     _S(declare_input_variable<Scalar>("effective_saturation")),
      56            2 :     _Pc(declare_output_variable<Scalar>("capillary_pressure")),
      57            4 :     _log_extension(options.get<bool>("log_extension")),
      58            4 :     _Sp_specified(options.user_specified("transition_saturation")),
      59            4 :     _Sp(options.get<double>("transition_saturation"))
      60              : {
      61            2 :   if (_log_extension && !_Sp_specified)
      62              :     throw NEMLException("BrooksCoreyCapillaryPressure: log_extension is set to true, but "
      63            0 :                         "transition_saturation is not specified.");
      64              : 
      65            2 :   if (!_log_extension && _Sp_specified)
      66              :     throw NEMLException("BrooksCoreyCapillaryPressure: transition_saturation is specified, but "
      67            0 :                         "log_extension is not set to true.");
      68              : 
      69            2 :   if (_Sp < 0.0 || _Sp > 1.0)
      70              :     throw NEMLException(
      71            0 :         "BrooksCoreyCapillaryPressure: transition_saturation must be in the range [0, 1].");
      72            2 : }
      73              : 
      74              : void
      75            6 : CapillaryPressure::set_value(bool out, bool dout, bool d2out)
      76              : {
      77            6 :   auto [Pc, dPc_dS, d2Pc_dS2] = calculate_pressure(_S, out, dout, d2out);
      78              : 
      79            6 :   if (_log_extension)
      80              :   {
      81              :     // Apply logarithmic extension
      82            6 :     constexpr double ln10 = 2.302585092994;
      83            6 :     auto [Pcs, dPcs_dS, d2Pcs_dS2] =
      84            6 :         calculate_pressure(Scalar::full(_Sp, _S.options()), true, true, false);
      85            6 :     auto slope = 1.0 / (ln10 * Pcs) * dPcs_dS;
      86            6 :     auto yintercept = log10(Pcs) - slope * _Sp;
      87            6 :     auto Pc_ext = Scalar(pow(10, slope * _S + yintercept));
      88              : 
      89            6 :     if (out)
      90            2 :       _Pc = where(_S < _Sp, Pc_ext, Pc);
      91              : 
      92            6 :     if (dout)
      93            2 :       _Pc.d(_S) = where(_S < _Sp, (ln10 * slope) * Pc_ext, dPc_dS);
      94              : 
      95            6 :     if (d2out)
      96            2 :       _Pc.d(_S, _S) = where(_S < _Sp, (ln10 * slope) * (ln10 * slope) * Pc_ext, d2Pc_dS2);
      97            6 :   }
      98              :   else
      99              :   {
     100            0 :     if (out)
     101            0 :       _Pc = Pc;
     102              : 
     103            0 :     if (dout)
     104            0 :       _Pc.d(_S) = dPc_dS;
     105              : 
     106            0 :     if (d2out)
     107            0 :       _Pc.d(_S, _S) = d2Pc_dS2;
     108              :   }
     109            6 : }
     110              : 
     111              : } // namespace neml2