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 "RichardsHalfGaussianSink.h"
      11             : 
      12             : // MOOSEincludes
      13             : #include "Function.h"
      14             : #include "Material.h"
      15             : #include "MooseVariable.h"
      16             : 
      17             : #include "libmesh/utility.h"
      18             : 
      19             : registerMooseObject("RichardsApp", RichardsHalfGaussianSink);
      20             : 
      21             : InputParameters
      22          22 : RichardsHalfGaussianSink::validParams()
      23             : {
      24          22 :   InputParameters params = IntegratedBC::validParams();
      25          44 :   params.addRequiredParam<Real>("max",
      26             :                                 "Maximum of the flux (measured in kg.m^-2.s^-1).  Flux out "
      27             :                                 "= max*exp((-0.5*(p - centre)/sd)^2) for p<centre, and Flux "
      28             :                                 "out = max for p>centre.  Note, to make this a source "
      29             :                                 "rather than a sink, let max<0");
      30          44 :   params.addRequiredParam<Real>("sd",
      31             :                                 "Standard deviation of the Gaussian (measured in Pa).  Flux "
      32             :                                 "out = max*exp((-0.5*(p - centre)/sd)^2) for p<centre, and "
      33             :                                 "Flux out = max for p>centre.");
      34          44 :   params.addRequiredParam<Real>("centre",
      35             :                                 "Centre of the Gaussian (measured in Pa).  Flux out = "
      36             :                                 "max*exp((-0.5*(p - centre)/sd)^2) for p<centre, and "
      37             :                                 "Flux out = max for p>centre.");
      38          44 :   params.addParam<FunctionName>("multiplying_fcn",
      39          44 :                                 1.0,
      40             :                                 "If this function is provided, the flux "
      41             :                                 "will be multiplied by this function.  "
      42             :                                 "This is useful for spatially or "
      43             :                                 "temporally varying sinks");
      44          44 :   params.addRequiredParam<UserObjectName>(
      45             :       "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
      46          22 :   return params;
      47           0 : }
      48             : 
      49          11 : RichardsHalfGaussianSink::RichardsHalfGaussianSink(const InputParameters & parameters)
      50             :   : IntegratedBC(parameters),
      51          11 :     _maximum(getParam<Real>("max")),
      52          22 :     _sd(getParam<Real>("sd")),
      53          22 :     _centre(getParam<Real>("centre")),
      54          11 :     _m_func(getFunction("multiplying_fcn")),
      55          11 :     _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
      56          11 :     _pvar(_richards_name_UO.richards_var_num(_var.number())),
      57          22 :     _pp(getMaterialProperty<std::vector<Real>>("porepressure")),
      58          33 :     _dpp_dv(getMaterialProperty<std::vector<std::vector<Real>>>("dporepressure_dv"))
      59             : {
      60          11 : }
      61             : 
      62             : Real
      63       75680 : RichardsHalfGaussianSink::computeQpResidual()
      64             : {
      65       75680 :   const Real test_fcn_f = _test[_i][_qp] * _m_func.value(_t, _q_point[_qp]);
      66             : 
      67       75680 :   if (_pp[_qp][_pvar] >= _centre)
      68        7520 :     return test_fcn_f * _maximum;
      69             : 
      70       68160 :   return test_fcn_f * _maximum *
      71       68160 :          std::exp(-0.5 * Utility::pow<2>((_pp[_qp][_pvar] - _centre) / _sd));
      72             : }
      73             : 
      74             : Real
      75      101632 : RichardsHalfGaussianSink::computeQpJacobian()
      76             : {
      77      101632 :   if (_pp[_qp][_pvar] >= _centre)
      78             :     return 0.0;
      79             : 
      80       81408 :   const Real test_fcn_f = _test[_i][_qp] * _m_func.value(_t, _q_point[_qp]);
      81       81408 :   return -test_fcn_f * _maximum * (_pp[_qp][_pvar] - _centre) / Utility::pow<2>(_sd) *
      82       81408 :          std::exp(-0.5 * Utility::pow<2>((_pp[_qp][_pvar] - _centre) / _sd)) * _phi[_j][_qp] *
      83       81408 :          _dpp_dv[_qp][_pvar][_pvar];
      84             : }
      85             : 
      86             : Real
      87       19200 : RichardsHalfGaussianSink::computeQpOffDiagJacobian(unsigned int jvar)
      88             : {
      89       19200 :   if (_richards_name_UO.not_richards_var(jvar))
      90             :     return 0.0;
      91             : 
      92       19200 :   if (_pp[_qp][_pvar] >= _centre)
      93             :     return 0.0;
      94             : 
      95       19200 :   const Real test_fcn_f = _test[_i][_qp] * _m_func.value(_t, _q_point[_qp]);
      96       19200 :   const unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
      97       19200 :   return -test_fcn_f * _maximum * (_pp[_qp][_pvar] - _centre) / Utility::pow<2>(_sd) *
      98       19200 :          std::exp(-0.5 * Utility::pow<2>((_pp[_qp][_pvar] - _centre) / _sd)) * _phi[_j][_qp] *
      99       19200 :          _dpp_dv[_qp][_pvar][dvar];
     100             : }