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             : //  This post processor returns the mass due to a half-gaussian sink flux from the boundary of a
      11             : //  volume.
      12             : //
      13             : #include "RichardsHalfGaussianSinkFlux.h"
      14             : #include "Function.h"
      15             : 
      16             : registerMooseObject("RichardsApp", RichardsHalfGaussianSinkFlux);
      17             : 
      18             : InputParameters
      19          16 : RichardsHalfGaussianSinkFlux::validParams()
      20             : {
      21          16 :   InputParameters params = SideIntegralVariablePostprocessor::validParams();
      22          32 :   params.addRequiredParam<Real>("max",
      23             :                                 "Maximum of the flux (measured in kg.m^-2.s^-1).  Flux out "
      24             :                                 "= max*exp((-0.5*(p - centre)/sd)^2) for p<centre, and Flux "
      25             :                                 "out = max for p>centre.  Note, to make this a source "
      26             :                                 "rather than a sink, let max<0");
      27          32 :   params.addRequiredParam<Real>("sd",
      28             :                                 "Standard deviation of the Gaussian (measured in Pa).  Flux "
      29             :                                 "out = max*exp((-0.5*(p - centre)/sd)^2) for p<centre, and "
      30             :                                 "Flux out = max for p>centre.");
      31          32 :   params.addRequiredParam<Real>("centre",
      32             :                                 "Centre of the Gaussian (measured in Pa).  Flux out = "
      33             :                                 "max*exp((-0.5*(p - centre)/sd)^2) for p<centre, and "
      34             :                                 "Flux out = max for p>centre.");
      35          32 :   params.addParam<FunctionName>(
      36             :       "multiplying_fcn",
      37          32 :       1.0,
      38             :       "The flux will be multiplied by this spatially-and-temporally varying function.");
      39          32 :   params.addRequiredParam<UserObjectName>(
      40             :       "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
      41          16 :   return params;
      42           0 : }
      43             : 
      44           8 : RichardsHalfGaussianSinkFlux::RichardsHalfGaussianSinkFlux(const InputParameters & parameters)
      45             :   : SideIntegralVariablePostprocessor(parameters),
      46           0 :     _feproblem(dynamic_cast<FEProblemBase &>(_subproblem)),
      47          16 :     _maximum(getParam<Real>("max")),
      48          16 :     _sd(getParam<Real>("sd")),
      49          16 :     _centre(getParam<Real>("centre")),
      50           8 :     _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
      51           8 :     _pvar(_richards_name_UO.richards_var_num(coupled("variable"))),
      52           8 :     _m_func(getFunction("multiplying_fcn")),
      53          24 :     _pp(getMaterialProperty<std::vector<Real>>("porepressure"))
      54             : {
      55           8 : }
      56             : 
      57             : Real
      58        1608 : RichardsHalfGaussianSinkFlux::computeQpIntegral()
      59             : {
      60        1608 :   if (_pp[_qp][_pvar] >= _centre)
      61         416 :     return _maximum * _dt * _m_func.value(_t, _q_point[_qp]);
      62             :   else
      63        1192 :     return _maximum * exp(-0.5 * std::pow((_pp[_qp][_pvar] - _centre) / _sd, 2)) * _dt *
      64        1192 :            _m_func.value(_t, _q_point[_qp]);
      65             : }