LCOV - code coverage report
Current view: top level - src/kernels - PorousFlowMassTimeDerivative.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 66 69 95.7 %
Date: 2017-11-20 14:50:56 Functions: 8 8 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /****************************************************************/
       2             : /* MOOSE - Multiphysics Object Oriented Simulation Environment  */
       3             : /*                                                              */
       4             : /*          All contents are licensed under LGPL V2.1           */
       5             : /*             See LICENSE for full restrictions                */
       6             : /****************************************************************/
       7             : 
       8             : #include "PorousFlowMassTimeDerivative.h"
       9             : 
      10             : // MOOSE includes
      11             : #include "MooseVariable.h"
      12             : 
      13             : #include "libmesh/quadrature.h"
      14             : 
      15             : // C++ includes
      16             : #include <limits>
      17             : 
      18             : template <>
      19             : InputParameters
      20         170 : validParams<PorousFlowMassTimeDerivative>()
      21             : {
      22         170 :   InputParameters params = validParams<TimeKernel>();
      23         510 :   params.addParam<bool>("strain_at_nearest_qp",
      24             :                         false,
      25             :                         "When calculating nodal porosity that depends on strain, use the strain at "
      26             :                         "the nearest quadpoint.  This adds a small extra computational burden, and "
      27             :                         "is not necessary for simulations involving only linear lagrange elements. "
      28             :                         " If you set this to true, you will also want to set the same parameter to "
      29         170 :                         "true for related Kernels and Materials");
      30         510 :   params.addParam<unsigned int>(
      31         170 :       "fluid_component", 0, "The index corresponding to the component for this kernel");
      32         510 :   params.addRequiredParam<UserObjectName>(
      33         170 :       "PorousFlowDictator", "The UserObject that holds the list of Porous-Flow variable names.");
      34         340 :   params.addClassDescription(
      35         170 :       "Component mass derivative wrt time for component given by fluid_component");
      36         170 :   return params;
      37             : }
      38             : 
      39         164 : PorousFlowMassTimeDerivative::PorousFlowMassTimeDerivative(const InputParameters & parameters)
      40             :   : TimeKernel(parameters),
      41         492 :     _fluid_component(getParam<unsigned int>("fluid_component")),
      42         328 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      43         164 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      44         164 :     _num_phases(_dictator.numPhases()),
      45         492 :     _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
      46         328 :     _porosity(getMaterialProperty<Real>("PorousFlow_porosity_nodal")),
      47         328 :     _porosity_old(getMaterialPropertyOld<Real>("PorousFlow_porosity_nodal")),
      48         328 :     _dporosity_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_nodal_dvar")),
      49             :     _dporosity_dgradvar(
      50         328 :         getMaterialProperty<std::vector<RealGradient>>("dPorousFlow_porosity_nodal_dgradvar")),
      51         164 :     _nearest_qp(_strain_at_nearest_qp
      52         329 :                     ? &getMaterialProperty<unsigned int>("PorousFlow_nearestqp_nodal")
      53             :                     : nullptr),
      54         328 :     _fluid_density(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
      55             :     _fluid_density_old(
      56         328 :         getMaterialPropertyOld<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
      57             :     _dfluid_density_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
      58         328 :         "dPorousFlow_fluid_phase_density_nodal_dvar")),
      59         328 :     _fluid_saturation_nodal(getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")),
      60             :     _fluid_saturation_nodal_old(
      61         328 :         getMaterialPropertyOld<std::vector<Real>>("PorousFlow_saturation_nodal")),
      62             :     _dfluid_saturation_nodal_dvar(
      63         328 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_saturation_nodal_dvar")),
      64         328 :     _mass_frac(getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
      65             :     _mass_frac_old(
      66         328 :         getMaterialPropertyOld<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
      67             :     _dmass_frac_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      68        3444 :         "dPorousFlow_mass_frac_nodal_dvar"))
      69             : {
      70         164 :   if (_fluid_component >= _dictator.numComponents())
      71           0 :     mooseError(
      72             :         "The Dictator proclaims that the number of components in this simulation is ",
      73           0 :         _dictator.numComponents(),
      74             :         " whereas you have used the Kernel PorousFlowComponetMassTimeDerivative with component = ",
      75             :         _fluid_component,
      76           0 :         ".  The Dictator does not take such mistakes lightly");
      77         164 : }
      78             : 
      79             : Real
      80     9547404 : PorousFlowMassTimeDerivative::computeQpResidual()
      81             : {
      82             :   Real mass = 0.0;
      83             :   Real mass_old = 0.0;
      84    31064772 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
      85             :   {
      86    43034736 :     mass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
      87    10758684 :             _mass_frac[_i][ph][_fluid_component];
      88    43034736 :     mass_old += _fluid_density_old[_i][ph] * _fluid_saturation_nodal_old[_i][ph] *
      89    10758684 :                 _mass_frac_old[_i][ph][_fluid_component];
      90             :   }
      91             : 
      92    38189616 :   return _test[_i][_qp] * (_porosity[_i] * mass - _porosity_old[_i] * mass_old) / _dt;
      93             : }
      94             : 
      95             : Real
      96    26777384 : PorousFlowMassTimeDerivative::computeQpJacobian()
      97             : {
      98             :   /// If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
      99    26777384 :   if (!_var_is_porflow_var)
     100             :     return 0.0;
     101    26777384 :   return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
     102             : }
     103             : 
     104             : Real
     105    21023072 : PorousFlowMassTimeDerivative::computeQpOffDiagJacobian(unsigned int jvar)
     106             : {
     107             :   /// If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
     108    21023072 :   if (_dictator.notPorousFlowVariable(jvar))
     109             :     return 0.0;
     110    21023072 :   return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     111             : }
     112             : 
     113             : Real
     114    47800456 : PorousFlowMassTimeDerivative::computeQpJac(unsigned int pvar)
     115             : {
     116    47806600 :   const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
     117             : 
     118             :   // porosity is dependent on variables that are lumped to the nodes,
     119             :   // but it can depend on the gradient
     120             :   // of variables, which are NOT lumped to the nodes, hence:
     121             :   Real dmass = 0.0;
     122   153309080 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     123   211017248 :     dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     124   158262936 :              _mass_frac[_i][ph][_fluid_component] * _dporosity_dgradvar[_i][pvar] *
     125    52754312 :              _grad_phi[_j][nearest_qp];
     126             : 
     127    47800456 :   if (_i != _j)
     128    79324088 :     return _test[_i][_qp] * dmass / _dt;
     129             : 
     130             :   /// As the fluid mass is lumped to the nodes, only non-zero terms are for _i==_j
     131    27456292 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     132             :   {
     133    48294700 :     dmass += _dfluid_density_dvar[_i][ph][pvar] * _fluid_saturation_nodal[_i][ph] *
     134    28976820 :              _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     135    38635760 :     dmass += _fluid_density[_i][ph] * _dfluid_saturation_nodal_dvar[_i][ph][pvar] *
     136     9658940 :              _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     137    28976820 :     dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     138    19317880 :              _dmass_frac_dvar[_i][ph][_fluid_component][pvar] * _porosity[_i];
     139    19317880 :     dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     140    19317880 :              _mass_frac[_i][ph][_fluid_component] * _dporosity_dvar[_i][pvar];
     141             :   }
     142    16276824 :   return _test[_i][_qp] * dmass / _dt;
     143        2499 : }

Generated by: LCOV version 1.11