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 "Q2PNodalMass.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseVariable.h"
      14             : 
      15             : // C++ includes
      16             : #include <iostream>
      17             : 
      18             : registerMooseObject("RichardsApp", Q2PNodalMass);
      19             : 
      20             : InputParameters
      21          70 : Q2PNodalMass::validParams()
      22             : {
      23          70 :   InputParameters params = TimeKernel::validParams();
      24         140 :   params.addRequiredParam<UserObjectName>(
      25             :       "fluid_density",
      26             :       "A RichardsDensity UserObject that defines the fluid density as a function of pressure.");
      27         140 :   params.addRequiredCoupledVar("other_var",
      28             :                                "The other variable in the 2-phase system.  If "
      29             :                                "Variable=porepressure, then other_var should be the "
      30             :                                "saturation Variable, and vice-versa.");
      31         140 :   params.addRequiredParam<bool>(
      32             :       "var_is_porepressure",
      33             :       "This flag is needed to correctly calculate the Jacobian entries.  If "
      34             :       "set to true, this Kernel will assume it is describing the mass of "
      35             :       "the phase with porepressure as its Variable (eg, the liquid phase).  "
      36             :       "If set to false, this Kernel will assumed it is describing the mass "
      37             :       "of the phase with saturation as its variable (eg, the gas phase)");
      38          70 :   params.addClassDescription("Fluid mass lumped to the nodes divided by dt");
      39          70 :   return params;
      40           0 : }
      41             : 
      42          28 : Q2PNodalMass::Q2PNodalMass(const InputParameters & parameters)
      43             :   : TimeKernel(parameters),
      44          28 :     _density(getUserObject<RichardsDensity>("fluid_density")),
      45          28 :     _other_var_nodal(coupledDofValues("other_var")),
      46          28 :     _other_var_num(coupled("other_var")),
      47          56 :     _var_is_pp(getParam<bool>("var_is_porepressure")),
      48          84 :     _porosity(getMaterialProperty<Real>("porosity"))
      49             : {
      50          28 : }
      51             : 
      52             : Real
      53       84384 : Q2PNodalMass::computeQpResidual()
      54             : {
      55             :   Real density;
      56             :   Real mass;
      57             : 
      58       84384 :   if (_var_is_pp)
      59             :   {
      60       42192 :     density = _density.density(_var.dofValues()[_i]);
      61       42192 :     mass = _porosity[_qp] * density * (1 - _other_var_nodal[_i]);
      62             :   }
      63             :   else
      64             :   {
      65       42192 :     density = _density.density(_other_var_nodal[_i]);
      66       42192 :     mass = _porosity[_qp] * density * _var.dofValues()[_i];
      67             :   }
      68             : 
      69       84384 :   return _test[_i][_qp] * mass / _dt;
      70             : }
      71             : 
      72             : Real
      73      213760 : Q2PNodalMass::computeQpJacobian()
      74             : {
      75      213760 :   if (_i != _j)
      76             :     return 0.0;
      77             : 
      78             :   Real mass_prime;
      79             : 
      80       46368 :   if (_var_is_pp)
      81             :   {
      82             :     // we're calculating the derivative wrt porepressure
      83       23184 :     Real ddensity = _density.ddensity(_var.dofValues()[_i]);
      84       23184 :     mass_prime = _porosity[_qp] * ddensity * (1 - _other_var_nodal[_i]);
      85             :   }
      86             :   else
      87             :   {
      88             :     // we're calculating the deriv wrt saturation
      89       23184 :     Real density = _density.density(_other_var_nodal[_i]);
      90       23184 :     mass_prime = _porosity[_qp] * density;
      91             :   }
      92             : 
      93       46368 :   return _test[_i][_qp] * mass_prime / _dt;
      94             : }
      95             : 
      96             : Real
      97      210560 : Q2PNodalMass::computeQpOffDiagJacobian(unsigned int jvar)
      98             : {
      99      210560 :   if (jvar != _other_var_num)
     100             :     return 0.0;
     101      209536 :   if (_i != _j)
     102             :     return 0.0;
     103             : 
     104             :   Real mass_prime;
     105             : 
     106       45312 :   if (_var_is_pp)
     107             :   {
     108             :     // we're calculating the deriv wrt saturation
     109       22656 :     Real density = _density.density(_var.dofValues()[_i]);
     110       22656 :     mass_prime = -_porosity[_qp] * density;
     111             :   }
     112             :   else
     113             :   {
     114             :     // we're calculating the deriv wrt porepressure
     115       22656 :     Real ddensity = _density.ddensity(_other_var_nodal[_i]);
     116       22656 :     mass_prime = _porosity[_qp] * ddensity * _var.dofValues()[_i];
     117             :   }
     118             : 
     119       45312 :   return _test[_i][_qp] * mass_prime / _dt;
     120             : }