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 "RichardsLumpedMassChange.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseVariable.h"
      14             : 
      15             : // C++ includes
      16             : #include <iostream>
      17             : 
      18             : registerMooseObject("RichardsApp", RichardsLumpedMassChange);
      19             : 
      20             : InputParameters
      21         122 : RichardsLumpedMassChange::validParams()
      22             : {
      23         122 :   InputParameters params = TimeKernel::validParams();
      24         244 :   params.addRequiredParam<UserObjectName>(
      25             :       "richardsVarNames_UO", "The UserObject that holds the list of Richards variables.");
      26         244 :   params.addRequiredParam<std::vector<UserObjectName>>(
      27             :       "density_UO",
      28             :       "List of names of user objects that define the fluid density (or densities for "
      29             :       "multiphase).  In the multiphase case, for ease of use, the density, Seff and "
      30             :       "Sat UserObjects are the same format as for RichardsMaterial, but only the one "
      31             :       "relevant for the specific phase is actually used.");
      32         244 :   params.addRequiredParam<std::vector<UserObjectName>>(
      33             :       "seff_UO",
      34             :       "List of name of user objects that define effective saturation as a function of "
      35             :       "porepressure(s)");
      36         244 :   params.addRequiredParam<std::vector<UserObjectName>>(
      37             :       "sat_UO",
      38             :       "List of names of user objects that define saturation as a function of effective saturation");
      39         122 :   return params;
      40           0 : }
      41             : 
      42          61 : RichardsLumpedMassChange::RichardsLumpedMassChange(const InputParameters & parameters)
      43             :   : TimeKernel(parameters),
      44          61 :     _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
      45          61 :     _num_p(_richards_name_UO.num_v()),
      46          61 :     _pvar(_richards_name_UO.richards_var_num(_var.number())),
      47             : 
      48         122 :     _porosity(getMaterialProperty<Real>("porosity")),
      49         122 :     _porosity_old(getMaterialProperty<Real>("porosity_old")),
      50             : 
      51             :     // in the following:  first get the userobject names that were inputted, then get the _pvar one
      52             :     // of these, then get the actual userobject that this corresponds to, then finally & gives
      53             :     // pointer to RichardsDensity (or whatever) object.
      54          61 :     _seff_UO(
      55         122 :         getUserObjectByName<RichardsSeff>(getParam<std::vector<UserObjectName>>("seff_UO")[_pvar])),
      56          61 :     _sat_UO(
      57         122 :         getUserObjectByName<RichardsSat>(getParam<std::vector<UserObjectName>>("sat_UO")[_pvar])),
      58          61 :     _density_UO(getUserObjectByName<RichardsDensity>(
      59         244 :         getParam<std::vector<UserObjectName>>("density_UO")[_pvar]))
      60             : {
      61          61 :   _ps_at_nodes.resize(_num_p);
      62          61 :   _ps_old_at_nodes.resize(_num_p);
      63             : 
      64         170 :   for (unsigned int pnum = 0; pnum < _num_p; ++pnum)
      65             :   {
      66         109 :     _ps_at_nodes[pnum] = _richards_name_UO.nodal_var(pnum);
      67         109 :     _ps_old_at_nodes[pnum] = _richards_name_UO.nodal_var_old(pnum);
      68             :   }
      69             : 
      70          61 :   _dseff.resize(_num_p);
      71          61 : }
      72             : 
      73             : Real
      74     5827832 : RichardsLumpedMassChange::computeQpResidual()
      75             : {
      76             :   // current values:
      77     5827832 :   const Real density = _density_UO.density((*_ps_at_nodes[_pvar])[_i]);
      78     5827832 :   const Real seff = _seff_UO.seff(_ps_at_nodes, _i);
      79     5827832 :   const Real sat = _sat_UO.sat(seff);
      80     5827832 :   const Real mass = _porosity[_qp] * density * sat;
      81             : 
      82             :   // old values:
      83     5827832 :   const Real density_old = _density_UO.density((*_ps_old_at_nodes[_pvar])[_i]);
      84     5827832 :   const Real seff_old = _seff_UO.seff(_ps_old_at_nodes, _i);
      85     5827832 :   const Real sat_old = _sat_UO.sat(seff_old);
      86     5827832 :   const Real mass_old = _porosity_old[_qp] * density_old * sat_old;
      87             : 
      88     5827832 :   return _test[_i][_qp] * (mass - mass_old) / _dt;
      89             : }
      90             : 
      91             : Real
      92    17244672 : RichardsLumpedMassChange::computeQpJacobian()
      93             : {
      94    17244672 :   if (_i != _j)
      95             :     return 0.0;
      96             : 
      97     2898944 :   const Real density = _density_UO.density((*_ps_at_nodes[_pvar])[_i]);
      98     2898944 :   const Real ddensity = _density_UO.ddensity((*_ps_at_nodes[_pvar])[_i]);
      99             : 
     100     2898944 :   const Real seff = _seff_UO.seff(_ps_at_nodes, _i);
     101     2898944 :   _seff_UO.dseff(_ps_at_nodes, _i, _dseff);
     102             : 
     103     2898944 :   const Real sat = _sat_UO.sat(seff);
     104     2898944 :   const Real dsat = _sat_UO.dsat(seff);
     105             : 
     106     2898944 :   const Real mass_prime = _porosity[_qp] * (ddensity * sat + density * _dseff[_pvar] * dsat);
     107             : 
     108     2898944 :   return _test[_i][_qp] * mass_prime / _dt;
     109             : }
     110             : 
     111             : Real
     112      537472 : RichardsLumpedMassChange::computeQpOffDiagJacobian(unsigned int jvar)
     113             : {
     114      537472 :   if (_richards_name_UO.not_richards_var(jvar))
     115             :     return 0.0;
     116      537472 :   if (_i != _j)
     117             :     return 0.0;
     118      267584 :   const unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
     119             : 
     120      267584 :   const Real density = _density_UO.density((*_ps_at_nodes[_pvar])[_i]);
     121             : 
     122      267584 :   const Real seff = _seff_UO.seff(_ps_at_nodes, _i);
     123      267584 :   _seff_UO.dseff(_ps_at_nodes, _i, _dseff);
     124             : 
     125      267584 :   const Real dsat = _sat_UO.dsat(seff);
     126             : 
     127      267584 :   const Real mass_prime = _porosity[_qp] * density * _dseff[dvar] * dsat;
     128             : 
     129      267584 :   return _test[_i][_qp] * mass_prime / _dt;
     130             : }