LCOV - code coverage report
Current view: top level - src/kernels - PorousFlowMassRadioactiveDecay.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 59 65 90.8 %
Date: 2017-11-17 17:48:31 Functions: 7 8 87.5 %
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 "PorousFlowMassRadioactiveDecay.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           1 : validParams<PorousFlowMassRadioactiveDecay>()
      21             : {
      22           1 :   InputParameters params = validParams<TimeKernel>();
      23           3 :   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           1 :                         "true for related Kernels and Materials");
      30           3 :   params.addParam<unsigned int>(
      31           1 :       "fluid_component", 0, "The index corresponding to the fluid component for this kernel");
      32           3 :   params.addRequiredParam<Real>("decay_rate",
      33           1 :                                 "The decay rate (units 1/time) for the fluid component");
      34           3 :   params.addRequiredParam<UserObjectName>(
      35           1 :       "PorousFlowDictator", "The UserObject that holds the list of Porous-Flow variable names.");
      36           2 :   params.addClassDescription("Radioactive decay of a fluid component");
      37           1 :   return params;
      38             : }
      39             : 
      40           1 : PorousFlowMassRadioactiveDecay::PorousFlowMassRadioactiveDecay(const InputParameters & parameters)
      41             :   : TimeKernel(parameters),
      42           3 :     _decay_rate(getParam<Real>("decay_rate")),
      43           3 :     _fluid_component(getParam<unsigned int>("fluid_component")),
      44           2 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      45           1 :     _var_is_porflow_var(_dictator.isPorousFlowVariable(_var.number())),
      46           1 :     _num_phases(_dictator.numPhases()),
      47           3 :     _strain_at_nearest_qp(getParam<bool>("strain_at_nearest_qp")),
      48           2 :     _porosity(getMaterialProperty<Real>("PorousFlow_porosity_nodal")),
      49           2 :     _dporosity_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_nodal_dvar")),
      50             :     _dporosity_dgradvar(
      51           2 :         getMaterialProperty<std::vector<RealGradient>>("dPorousFlow_porosity_nodal_dgradvar")),
      52           1 :     _nearest_qp(_strain_at_nearest_qp
      53           2 :                     ? &getMaterialProperty<unsigned int>("PorousFlow_nearestqp_nodal")
      54             :                     : nullptr),
      55           2 :     _fluid_density(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")),
      56             :     _dfluid_density_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
      57           2 :         "dPorousFlow_fluid_phase_density_nodal_dvar")),
      58           2 :     _fluid_saturation_nodal(getMaterialProperty<std::vector<Real>>("PorousFlow_saturation_nodal")),
      59             :     _dfluid_saturation_nodal_dvar(
      60           2 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_saturation_nodal_dvar")),
      61           2 :     _mass_frac(getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")),
      62             :     _dmass_frac_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      63          18 :         "dPorousFlow_mass_frac_nodal_dvar"))
      64             : {
      65           1 :   if (_fluid_component >= _dictator.numComponents())
      66           0 :     mooseError("The Dictator proclaims that the number of components in this simulation is ",
      67           0 :                _dictator.numComponents(),
      68             :                " whereas you have used the Kernel PorousFlowComponetMassRadioactiveDecay with "
      69             :                "component = ",
      70             :                _fluid_component,
      71           0 :                ".  The Dictator does not take such mistakes lightly");
      72           1 : }
      73             : 
      74             : Real
      75          84 : PorousFlowMassRadioactiveDecay::computeQpResidual()
      76             : {
      77             :   Real mass = 0.0;
      78         252 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
      79         336 :     mass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
      80          84 :             _mass_frac[_i][ph][_fluid_component];
      81             : 
      82         252 :   return _test[_i][_qp] * _decay_rate * _porosity[_i] * mass;
      83             : }
      84             : 
      85             : Real
      86         120 : PorousFlowMassRadioactiveDecay::computeQpJacobian()
      87             : {
      88             :   /// If the variable is not a PorousFlow variable (very unusual), the diag Jacobian terms are 0
      89         120 :   if (!_var_is_porflow_var)
      90             :     return 0.0;
      91         120 :   return computeQpJac(_dictator.porousFlowVariableNum(_var.number()));
      92             : }
      93             : 
      94             : Real
      95           0 : PorousFlowMassRadioactiveDecay::computeQpOffDiagJacobian(unsigned int jvar)
      96             : {
      97             :   /// If the variable is not a PorousFlow variable, the OffDiag Jacobian terms are 0
      98           0 :   if (_dictator.notPorousFlowVariable(jvar))
      99             :     return 0.0;
     100           0 :   return computeQpJac(_dictator.porousFlowVariableNum(jvar));
     101             : }
     102             : 
     103             : Real
     104         120 : PorousFlowMassRadioactiveDecay::computeQpJac(unsigned int pvar)
     105             : {
     106         120 :   const unsigned nearest_qp = (_strain_at_nearest_qp ? (*_nearest_qp)[_i] : _i);
     107             : 
     108             :   // porosity is dependent on variables that are lumped to the nodes,
     109             :   // but it can depend on the gradient
     110             :   // of variables, which are NOT lumped to the nodes, hence:
     111             :   Real dmass = 0.0;
     112         360 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     113         480 :     dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     114         360 :              _mass_frac[_i][ph][_fluid_component] * _dporosity_dgradvar[_i][pvar] *
     115         120 :              _grad_phi[_j][nearest_qp];
     116             : 
     117         120 :   if (_i != _j)
     118         120 :     return _test[_i][_qp] * _decay_rate * dmass;
     119             : 
     120             :   /// As the fluid mass is lumped to the nodes, only non-zero terms are for _i==_j
     121         180 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
     122             :   {
     123         300 :     dmass += _dfluid_density_dvar[_i][ph][pvar] * _fluid_saturation_nodal[_i][ph] *
     124         180 :              _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     125         240 :     dmass += _fluid_density[_i][ph] * _dfluid_saturation_nodal_dvar[_i][ph][pvar] *
     126          60 :              _mass_frac[_i][ph][_fluid_component] * _porosity[_i];
     127         180 :     dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     128         120 :              _dmass_frac_dvar[_i][ph][_fluid_component][pvar] * _porosity[_i];
     129         120 :     dmass += _fluid_density[_i][ph] * _fluid_saturation_nodal[_i][ph] *
     130         120 :              _mass_frac[_i][ph][_fluid_component] * _dporosity_dvar[_i][pvar];
     131             :   }
     132         120 :   return _test[_i][_qp] * _decay_rate * dmass;
     133        2499 : }

Generated by: LCOV version 1.11