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 "RichardsPPenalty.h"
      11             : 
      12             : #include <iostream>
      13             : 
      14             : registerMooseObject("RichardsApp", RichardsPPenalty);
      15             : 
      16             : InputParameters
      17          24 : RichardsPPenalty::validParams()
      18             : {
      19          24 :   InputParameters params = Kernel::validParams();
      20          48 :   params.addParam<Real>(
      21             :       "a",
      22          48 :       1.0E-10,
      23             :       "Weight of the penalty.  Penalty = a*(lower - variable) for variable<lower, "
      24             :       "and zero otherwise.  Care should be taken with this parameter choice.  "
      25             :       "Determine the typical size of your residual (usually rho*perm*(gradP - "
      26             :       "rho*g)/visc), then typically you want the penalty to ensure p>lower*(1-1E-6), "
      27             :       "so for the PPP formulation you typically Penalty = a*1E-6*|p|.  I recommend "
      28             :       "that Penalty = 1E-3*residual, yielding a = 1E3*residual/|P|. ");
      29          48 :   params.addRequiredCoupledVar(
      30             :       "lower_var", "Your variable will be constrained to be greater than this lower_var variable.");
      31          24 :   params.addClassDescription("This adds a term to the residual that attempts to enforce variable > "
      32             :                              "lower_var.  The term is a*(lower - variable) for variable<lower, and "
      33             :                              "zero otherwise");
      34          24 :   return params;
      35           0 : }
      36             : 
      37          12 : RichardsPPenalty::RichardsPPenalty(const InputParameters & parameters)
      38             :   : Kernel(parameters),
      39          12 :     _a(getParam<Real>("a")),
      40          12 :     _lower(coupledValue("lower_var")),
      41          24 :     _lower_var_num(coupled("lower_var"))
      42             : {
      43          12 : }
      44             : 
      45             : Real
      46       29000 : RichardsPPenalty::computeQpResidual()
      47             : {
      48       29000 :   if (_u[_qp] < _lower[_qp])
      49       17608 :     return _test[_i][_qp] * _a * (_lower[_qp] - _u[_qp]);
      50             : 
      51             :   return 0.0;
      52             : }
      53             : 
      54             : Real
      55       24320 : RichardsPPenalty::computeQpJacobian()
      56             : {
      57       24320 :   if (_u[_qp] < _lower[_qp])
      58       10764 :     return -_test[_i][_qp] * _a * _phi[_j][_qp];
      59             :   ;
      60             : 
      61             :   return 0.0;
      62             : }
      63             : 
      64             : Real
      65       24320 : RichardsPPenalty::computeQpOffDiagJacobian(unsigned int jvar)
      66             : {
      67       24320 :   if (jvar == _lower_var_num && _u[_qp] < _lower[_qp])
      68       10764 :     return _test[_i][_qp] * _a * _phi[_j][_qp];
      69             : 
      70             :   return 0.0;
      71             : }