LCOV - code coverage report
Current view: top level - src/kernels - PorousFlowDispersiveFlux.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 114 118 96.6 %
Date: 2017-11-21 14:47:27 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 "PorousFlowDispersiveFlux.h"
       9             : 
      10             : // MOOSE includes
      11             : #include "MooseVariable.h"
      12             : 
      13             : template <>
      14             : InputParameters
      15          20 : validParams<PorousFlowDispersiveFlux>()
      16             : {
      17          20 :   InputParameters params = validParams<Kernel>();
      18          60 :   params.addParam<unsigned int>(
      19          20 :       "fluid_component", 0, "The index corresponding to the fluid component for this kernel");
      20          60 :   params.addRequiredParam<UserObjectName>(
      21          20 :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      22          60 :   params.addRequiredParam<std::vector<Real>>(
      23          20 :       "disp_long", "Vector of longitudinal dispersion coefficients for each phase");
      24          60 :   params.addRequiredParam<std::vector<Real>>(
      25          20 :       "disp_trans", "Vector of transverse dispersion coefficients for each phase");
      26          60 :   params.addRequiredParam<RealVectorValue>("gravity",
      27          20 :                                            "Gravitational acceleration vector downwards (m/s^2)");
      28          40 :   params.addClassDescription(
      29          20 :       "Dispersive and diffusive flux of the component given by fluid_component in all phases");
      30          20 :   return params;
      31             : }
      32             : 
      33          20 : PorousFlowDispersiveFlux::PorousFlowDispersiveFlux(const InputParameters & parameters)
      34             :   : Kernel(parameters),
      35             : 
      36          40 :     _fluid_density_qp(getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_qp")),
      37             :     _dfluid_density_qp_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
      38          40 :         "dPorousFlow_fluid_phase_density_qp_dvar")),
      39             :     _grad_mass_frac(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
      40          40 :         "PorousFlow_grad_mass_frac_qp")),
      41             :     _dmass_frac_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      42          40 :         "dPorousFlow_mass_frac_qp_dvar")),
      43          40 :     _porosity_qp(getMaterialProperty<Real>("PorousFlow_porosity_qp")),
      44          40 :     _dporosity_qp_dvar(getMaterialProperty<std::vector<Real>>("dPorousFlow_porosity_qp_dvar")),
      45          40 :     _tortuosity(getMaterialProperty<std::vector<Real>>("PorousFlow_tortuosity_qp")),
      46             :     _dtortuosity_dvar(
      47          40 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_tortuosity_qp_dvar")),
      48             :     _diffusion_coeff(
      49          40 :         getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_diffusion_coeff_qp")),
      50             :     _ddiffusion_coeff_dvar(getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      51          40 :         "dPorousFlow_diffusion_coeff_qp_dvar")),
      52          40 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      53          60 :     _fluid_component(getParam<unsigned int>("fluid_component")),
      54          20 :     _num_phases(_dictator.numPhases()),
      55             :     _identity_tensor(RankTwoTensor::initIdentity),
      56             :     _relative_permeability(
      57          40 :         getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_qp")),
      58             :     _drelative_permeability_dvar(getMaterialProperty<std::vector<std::vector<Real>>>(
      59          40 :         "dPorousFlow_relative_permeability_qp_dvar")),
      60          40 :     _fluid_viscosity(getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_qp")),
      61             :     _dfluid_viscosity_dvar(
      62          40 :         getMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_qp_dvar")),
      63          40 :     _permeability(getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")),
      64             :     _dpermeability_dvar(
      65          40 :         getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
      66             :     _dpermeability_dgradvar(getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
      67          40 :         "dPorousFlow_permeability_qp_dgradvar")),
      68          40 :     _grad_p(getMaterialProperty<std::vector<RealGradient>>("PorousFlow_grad_porepressure_qp")),
      69             :     _dgrad_p_dgrad_var(getMaterialProperty<std::vector<std::vector<Real>>>(
      70          40 :         "dPorousFlow_grad_porepressure_qp_dgradvar")),
      71             :     _dgrad_p_dvar(getMaterialProperty<std::vector<std::vector<RealGradient>>>(
      72          40 :         "dPorousFlow_grad_porepressure_qp_dvar")),
      73          40 :     _gravity(getParam<RealVectorValue>("gravity")),
      74             :     _disp_long(getParam<std::vector<Real>>("disp_long")),
      75         600 :     _disp_trans(getParam<std::vector<Real>>("disp_trans"))
      76             : {
      77             :   // Check that sufficient values of the dispersion coefficients have been entered
      78          20 :   if (_disp_long.size() != _num_phases)
      79             :     mooseError("The number of longitudinal dispersion coefficients disp_long in ",
      80           0 :                _name,
      81           0 :                " is not equal to the number of phases");
      82             : 
      83          20 :   if (_disp_trans.size() != _num_phases)
      84             :     mooseError("The number of transverse dispersion coefficients disp_trans in ",
      85           0 :                _name,
      86           0 :                " is not equal to the number of phases");
      87          20 : }
      88             : 
      89             : Real
      90      125120 : PorousFlowDispersiveFlux::computeQpResidual()
      91             : {
      92             :   RealVectorValue flux = 0.0;
      93             :   RealVectorValue velocity;
      94             :   Real velocity_abs;
      95      125120 :   RankTwoTensor v2;
      96      125120 :   RankTwoTensor dispersion;
      97      125120 :   dispersion.zero();
      98             :   Real diffusion;
      99             : 
     100      396480 :   for (unsigned int ph = 0; ph < _num_phases; ++ph)
     101             :   {
     102             :     // Diffusive component
     103      135680 :     diffusion =
     104      542720 :         _porosity_qp[_qp] * _tortuosity[_qp][ph] * _diffusion_coeff[_qp][ph][_fluid_component];
     105             : 
     106             :     // Calculate Darcy velocity
     107      678400 :     velocity = (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity) *
     108      271360 :                 _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph]);
     109      135680 :     velocity_abs = std::sqrt(velocity * velocity);
     110             : 
     111      135680 :     if (velocity_abs > 0.0)
     112             :     {
     113      117384 :       v2.vectorOuterProduct(velocity, velocity);
     114             : 
     115             :       // Add longitudinal dispersion to diffusive component
     116      117384 :       diffusion += _disp_trans[ph] * velocity_abs;
     117      234768 :       dispersion = (_disp_long[ph] - _disp_trans[ph]) * v2 / velocity_abs;
     118             :     }
     119             : 
     120      678400 :     flux += _fluid_density_qp[_qp][ph] * (diffusion * _identity_tensor + dispersion) *
     121      135680 :             _grad_mass_frac[_qp][ph][_fluid_component];
     122             :   }
     123      250240 :   return _grad_test[_i][_qp] * flux;
     124             : }
     125             : 
     126             : Real
     127       69184 : PorousFlowDispersiveFlux::computeQpJacobian()
     128             : {
     129       69184 :   return computeQpJac(_var.number());
     130             : }
     131             : 
     132             : Real
     133       69184 : PorousFlowDispersiveFlux::computeQpOffDiagJacobian(unsigned int jvar)
     134             : {
     135       69184 :   return computeQpJac(jvar);
     136             : }
     137             : 
     138             : Real
     139      138368 : PorousFlowDispersiveFlux::computeQpJac(unsigned int jvar) const
     140             : {
     141             :   // If the variable is not a valid PorousFlow variable, set the Jacobian to 0
     142      138368 :   if (_dictator.notPorousFlowVariable(jvar))
     143             :     return 0.0;
     144             : 
     145      138368 :   const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
     146             : 
     147             :   RealVectorValue velocity;
     148             :   Real velocity_abs;
     149      138368 :   RankTwoTensor v2;
     150      138368 :   RankTwoTensor dispersion;
     151      138368 :   dispersion.zero();
     152             :   Real diffusion;
     153             :   RealVectorValue flux = 0.0;
     154             :   RealVectorValue dflux = 0.0;
     155             : 
     156      424320 :   for (unsigned int ph = 0; ph < _num_phases; ++ph)
     157             :   {
     158             :     // Diffusive component
     159      142976 :     diffusion =
     160      571904 :         _porosity_qp[_qp] * _tortuosity[_qp][ph] * _diffusion_coeff[_qp][ph][_fluid_component];
     161             : 
     162             :     // Calculate Darcy velocity
     163      714880 :     velocity = (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity) *
     164      285952 :                 _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph]);
     165      142976 :     velocity_abs = std::sqrt(velocity * velocity);
     166             : 
     167      142976 :     if (velocity_abs > 0.0)
     168             :     {
     169      122656 :       v2.vectorOuterProduct(velocity, velocity);
     170             : 
     171             :       // Add longitudinal dispersion to diffusive component
     172      122656 :       diffusion += _disp_trans[ph] * velocity_abs;
     173      245312 :       dispersion = (_disp_long[ph] - _disp_trans[ph]) * v2 / velocity_abs;
     174             :     }
     175             : 
     176             :     // Derivative of Darcy velocity
     177      428928 :     RealVectorValue dvelocity = _dpermeability_dvar[_qp][pvar] * _phi[_j][_qp] *
     178      571904 :                                 (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
     179     1000832 :     for (unsigned i = 0; i < LIBMESH_DIM; ++i)
     180     1286784 :       dvelocity += _dpermeability_dgradvar[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
     181      857856 :                    (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
     182             :     dvelocity +=
     183      571904 :         _permeability[_qp] * (_grad_phi[_j][_qp] * _dgrad_p_dgrad_var[_qp][ph][pvar] -
     184      285952 :                               _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] * _gravity);
     185      285952 :     dvelocity += _permeability[_qp] * (_dgrad_p_dvar[_qp][ph][pvar] * _phi[_j][_qp]);
     186             : 
     187             :     Real dvelocity_abs = 0.0;
     188      142976 :     if (velocity_abs > 0.0)
     189      122656 :       dvelocity_abs = velocity * dvelocity / velocity_abs;
     190             : 
     191             :     // Derivative of diffusion term (note: dispersivity is assumed constant)
     192      571904 :     Real ddiffusion = _phi[_j][_qp] * _dporosity_qp_dvar[_qp][pvar] * _tortuosity[_qp][ph] *
     193      285952 :                       _diffusion_coeff[_qp][ph][_fluid_component];
     194      428928 :     ddiffusion += _phi[_j][_qp] * _porosity_qp[_qp] * _dtortuosity_dvar[_qp][ph][pvar] *
     195             :                   _diffusion_coeff[_qp][ph][_fluid_component];
     196      285952 :     ddiffusion += _phi[_j][_qp] * _porosity_qp[_qp] * _tortuosity[_qp][ph] *
     197      142976 :                   _ddiffusion_coeff_dvar[_qp][ph][_fluid_component][pvar];
     198      142976 :     ddiffusion += _disp_trans[ph] * dvelocity_abs;
     199             : 
     200             :     // Derivative of dispersion term (note: dispersivity is assumed constant)
     201      142976 :     RankTwoTensor ddispersion;
     202      142976 :     ddispersion.zero();
     203      142976 :     if (velocity_abs > 0.0)
     204             :     {
     205      122656 :       RankTwoTensor dv2a, dv2b;
     206      122656 :       dv2a.vectorOuterProduct(velocity, dvelocity);
     207      122656 :       dv2b.vectorOuterProduct(dvelocity, velocity);
     208      490624 :       ddispersion = (_disp_long[ph] - _disp_trans[ph]) * (dv2a + dv2b) / velocity_abs;
     209             :       ddispersion -=
     210      367968 :           (_disp_long[ph] - _disp_trans[ph]) * v2 * dvelocity_abs / velocity_abs / velocity_abs;
     211             :     }
     212             : 
     213      428928 :     dflux += _phi[_j][_qp] * _dfluid_density_qp_dvar[_qp][ph][pvar] *
     214     1000832 :              (diffusion * _identity_tensor + dispersion) *
     215      142976 :              _grad_mass_frac[_qp][ph][_fluid_component];
     216      714880 :     dflux += _fluid_density_qp[_qp][ph] * (ddiffusion * _identity_tensor + ddispersion) *
     217      142976 :              _grad_mass_frac[_qp][ph][_fluid_component];
     218      571904 :     dflux += _fluid_density_qp[_qp][ph] * (diffusion * _identity_tensor + dispersion) *
     219      714880 :              _dmass_frac_dvar[_qp][ph][_fluid_component][pvar] * _grad_phi[_j][_qp];
     220             :   }
     221             : 
     222      138368 :   return _grad_test[_i][_qp] * dflux;
     223        2499 : }

Generated by: LCOV version 1.11