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 "PorousFlowDispersiveFlux.h"
      11             : 
      12             : #include "MooseVariable.h"
      13             : 
      14             : registerMooseObject("PorousFlowApp", PorousFlowDispersiveFlux);
      15             : registerMooseObject("PorousFlowApp", ADPorousFlowDispersiveFlux);
      16             : 
      17             : template <bool is_ad>
      18             : InputParameters
      19        1376 : PorousFlowDispersiveFluxTempl<is_ad>::validParams()
      20             : {
      21             :   InputParameters params = GenericKernel<is_ad>::validParams();
      22        1402 :   params.addParam<unsigned int>(
      23        1402 :       "fluid_component", 0, "The index corresponding to the fluid component for this kernel");
      24        1402 :   params.addRequiredParam<UserObjectName>(
      25             :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      26        1402 :   params.addRequiredParam<std::vector<Real>>(
      27             :       "disp_long", "Vector of longitudinal dispersion coefficients for each phase");
      28        1402 :   params.addRequiredParam<std::vector<Real>>(
      29             :       "disp_trans", "Vector of transverse dispersion coefficients for each phase");
      30        1402 :   params.addRequiredParam<RealVectorValue>("gravity",
      31             :                                            "Gravitational acceleration vector downwards (m/s^2)");
      32         701 :   params.addClassDescription(
      33             :       "Dispersive and diffusive flux of the component given by fluid_component in all phases");
      34         701 :   return params;
      35           0 : }
      36             : 
      37             : template <bool is_ad>
      38         372 : PorousFlowDispersiveFluxTempl<is_ad>::PorousFlowDispersiveFluxTempl(
      39             :     const InputParameters & parameters)
      40             :   : GenericKernel<is_ad>(parameters),
      41         372 :     _fluid_density_qp(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      42             :         "PorousFlow_fluid_phase_density_qp")),
      43         372 :     _dfluid_density_qp_dvar(
      44             :         is_ad ? nullptr
      45         358 :               : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      46             :                     "dPorousFlow_fluid_phase_density_qp_dvar")),
      47         372 :     _grad_mass_frac(
      48         372 :         this->template getGenericMaterialProperty<std::vector<std::vector<RealGradient>>, is_ad>(
      49             :             "PorousFlow_grad_mass_frac_qp")),
      50         372 :     _dmass_frac_dvar(
      51             :         is_ad ? nullptr
      52         358 :               : &this->template getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      53             :                     "dPorousFlow_mass_frac_qp_dvar")),
      54         744 :     _porosity_qp(this->template getGenericMaterialProperty<Real, is_ad>("PorousFlow_porosity_qp")),
      55         372 :     _dporosity_qp_dvar(is_ad ? nullptr
      56         358 :                              : &this->template getMaterialProperty<std::vector<Real>>(
      57             :                                    "dPorousFlow_porosity_qp_dvar")),
      58         744 :     _tortuosity(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      59             :         "PorousFlow_tortuosity_qp")),
      60         372 :     _dtortuosity_dvar(is_ad ? nullptr
      61         358 :                             : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      62             :                                   "dPorousFlow_tortuosity_qp_dvar")),
      63         744 :     _diffusion_coeff(this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      64             :         "PorousFlow_diffusion_coeff_qp")),
      65         372 :     _ddiffusion_coeff_dvar(
      66             :         is_ad ? nullptr
      67         358 :               : &this->template getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      68             :                     "dPorousFlow_diffusion_coeff_qp_dvar")),
      69         372 :     _dictator(this->template getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      70         744 :     _fluid_component(this->template getParam<unsigned int>("fluid_component")),
      71         372 :     _num_phases(_dictator.numPhases()),
      72         372 :     _identity_tensor(GenericRankTwoTensor<is_ad>::initIdentity),
      73         744 :     _relative_permeability(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      74             :         "PorousFlow_relative_permeability_qp")),
      75         372 :     _drelative_permeability_dvar(
      76             :         is_ad ? nullptr
      77         358 :               : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      78             :                     "dPorousFlow_relative_permeability_qp_dvar")),
      79         744 :     _fluid_viscosity(this->template getGenericMaterialProperty<std::vector<Real>, is_ad>(
      80             :         "PorousFlow_viscosity_qp")),
      81         372 :     _dfluid_viscosity_dvar(
      82             :         is_ad ? nullptr
      83         358 :               : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      84             :                     "dPorousFlow_viscosity_qp_dvar")),
      85         744 :     _permeability(this->template getGenericMaterialProperty<RealTensorValue, is_ad>(
      86             :         "PorousFlow_permeability_qp")),
      87         372 :     _dpermeability_dvar(is_ad ? nullptr
      88         358 :                               : &this->template getMaterialProperty<std::vector<RealTensorValue>>(
      89             :                                     "dPorousFlow_permeability_qp_dvar")),
      90         372 :     _dpermeability_dgradvar(
      91             :         is_ad ? nullptr
      92         358 :               : &this->template getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
      93             :                     "dPorousFlow_permeability_qp_dgradvar")),
      94         744 :     _grad_p(this->template getGenericMaterialProperty<std::vector<RealGradient>, is_ad>(
      95             :         "PorousFlow_grad_porepressure_qp")),
      96         372 :     _dgrad_p_dgrad_var(is_ad ? nullptr
      97         358 :                              : &this->template getMaterialProperty<std::vector<std::vector<Real>>>(
      98             :                                    "dPorousFlow_grad_porepressure_qp_dgradvar")),
      99         372 :     _dgrad_p_dvar(is_ad
     100             :                       ? nullptr
     101         358 :                       : &this->template getMaterialProperty<std::vector<std::vector<RealGradient>>>(
     102             :                             "dPorousFlow_grad_porepressure_qp_dvar")),
     103         744 :     _gravity(this->template getParam<RealVectorValue>("gravity")),
     104         744 :     _disp_long(this->template getParam<std::vector<Real>>("disp_long")),
     105         744 :     _disp_trans(this->template getParam<std::vector<Real>>("disp_trans")),
     106         744 :     _perm_derivs(_dictator.usePermDerivs())
     107             : {
     108         372 :   if (_fluid_component >= _dictator.numComponents())
     109           0 :     this->paramError(
     110             :         "fluid_component",
     111             :         "The Dictator proclaims that the maximum fluid component index in this simulation is ",
     112           0 :         _dictator.numComponents() - 1,
     113             :         " whereas you have used ",
     114           0 :         _fluid_component,
     115             :         ". Remember that indexing starts at 0. The Dictator does not take such mistakes lightly.");
     116             : 
     117         372 :   if (_disp_long.size() != _num_phases)
     118           0 :     this->paramError(
     119             :         "disp_long",
     120             :         "The number of longitudinal dispersion coefficients is not equal to the number of phases");
     121             : 
     122         372 :   if (_disp_trans.size() != _num_phases)
     123           0 :     this->paramError("disp_trans",
     124             :                      "The number of transverse dispersion coefficients disp_trans is not equal to "
     125             :                      "the number of phases");
     126         372 : }
     127             : 
     128             : template <bool is_ad>
     129             : GenericReal<is_ad>
     130     3767344 : PorousFlowDispersiveFluxTempl<is_ad>::computeQpResidual()
     131             : {
     132             :   GenericRealVectorValue<is_ad> flux;
     133             :   GenericRealVectorValue<is_ad> velocity;
     134             :   GenericReal<is_ad> velocity_abs;
     135     3767344 :   GenericRankTwoTensor<is_ad> v2;
     136     3767344 :   GenericRankTwoTensor<is_ad> dispersion;
     137        8064 :   dispersion.zero();
     138             :   GenericReal<is_ad> diffusion;
     139             : 
     140     7708688 :   for (unsigned int ph = 0; ph < _num_phases; ++ph)
     141             :   {
     142     3941344 :     diffusion =
     143     3941344 :         _porosity_qp[_qp] * _tortuosity[_qp][ph] * _diffusion_coeff[_qp][ph][_fluid_component];
     144             : 
     145     3941344 :     velocity = _permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity) *
     146     3941344 :                _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph];
     147     3941344 :     velocity_abs = velocity.norm();
     148             : 
     149     3941344 :     if (MetaPhysicL::raw_value(velocity_abs) > 0.0)
     150             :     {
     151     3478136 :       v2 = GenericRankTwoTensor<is_ad>::selfOuterProduct(velocity);
     152     3478136 :       diffusion += GenericReal<is_ad>(_disp_trans[ph]) * velocity_abs;
     153     3486200 :       dispersion = GenericReal<is_ad>(_disp_long[ph] - _disp_trans[ph]) * v2 / velocity_abs;
     154             :     }
     155             : 
     156     3949408 :     flux += _fluid_density_qp[_qp][ph] * (diffusion * _identity_tensor + dispersion) *
     157     3941344 :             _grad_mass_frac[_qp][ph][_fluid_component];
     158             :   }
     159     3767344 :   return _grad_test[_i][_qp] * flux;
     160             : }
     161             : 
     162             : template <bool is_ad>
     163             : Real
     164     9614832 : PorousFlowDispersiveFluxTempl<is_ad>::computeQpJacobian()
     165             : {
     166     9614832 :   return computeQpJac(_var.number());
     167             : }
     168             : 
     169             : template <bool is_ad>
     170             : Real
     171     9614272 : PorousFlowDispersiveFluxTempl<is_ad>::computeQpOffDiagJacobian(unsigned int jvar)
     172             : {
     173     9614272 :   return computeQpJac(jvar);
     174             : }
     175             : 
     176             : template <bool is_ad>
     177             : Real
     178    19229104 : PorousFlowDispersiveFluxTempl<is_ad>::computeQpJac(unsigned int jvar) const
     179             : {
     180             :   if constexpr (!is_ad)
     181             :   {
     182    19229104 :     if (_dictator.notPorousFlowVariable(jvar))
     183             :       return 0.0;
     184             : 
     185    19229104 :     const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
     186             : 
     187             :     RealVectorValue velocity;
     188             :     Real velocity_abs;
     189    19229104 :     RankTwoTensor v2;
     190    19229104 :     RankTwoTensor dispersion;
     191             :     dispersion.zero();
     192             :     Real diffusion;
     193             :     RealVectorValue dflux = 0.0;
     194             : 
     195    39573392 :     for (unsigned int ph = 0; ph < _num_phases; ++ph)
     196             :     {
     197    20344288 :       diffusion =
     198    20344288 :           _porosity_qp[_qp] * _tortuosity[_qp][ph] * _diffusion_coeff[_qp][ph][_fluid_component];
     199             : 
     200    20344288 :       velocity = _permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity) *
     201    20344288 :                  _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph];
     202    20344288 :       velocity_abs = MetaPhysicL::raw_value(velocity).norm();
     203             : 
     204    20344288 :       if (velocity_abs > 0.0)
     205             :       {
     206    18202624 :         v2 = RankTwoTensor::selfOuterProduct(velocity);
     207    18202624 :         diffusion += _disp_trans[ph] * velocity_abs;
     208    18202624 :         dispersion = (_disp_long[ph] - _disp_trans[ph]) * v2 / velocity_abs;
     209             :       }
     210             : 
     211             :       RealVectorValue dvelocity =
     212    20344288 :           _permeability[_qp] *
     213    20344288 :           (_grad_phi[_j][_qp] * (*_dgrad_p_dgrad_var)[_qp][ph][pvar] -
     214    20344288 :            _phi[_j][_qp] * (*_dfluid_density_qp_dvar)[_qp][ph][pvar] * _gravity);
     215    20344288 :       dvelocity += _permeability[_qp] * ((*_dgrad_p_dvar)[_qp][ph][pvar] * _phi[_j][_qp]);
     216             : 
     217    20344288 :       if (_perm_derivs)
     218             :       {
     219           0 :         dvelocity += (*_dpermeability_dvar)[_qp][pvar] * _phi[_j][_qp] *
     220           0 :                      (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
     221             : 
     222           0 :         for (const auto i : make_range(Moose::dim))
     223           0 :           dvelocity += (*_dpermeability_dgradvar)[_qp][i][pvar] * _grad_phi[_j][_qp](i) *
     224           0 :                        (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity);
     225             :       }
     226             : 
     227    20344288 :       dvelocity =
     228             :           dvelocity * _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph] +
     229    20344288 :           (_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density_qp[_qp][ph] * _gravity)) *
     230    20344288 :               ((*_drelative_permeability_dvar)[_qp][ph][pvar] / _fluid_viscosity[_qp][ph] -
     231    20344288 :                _relative_permeability[_qp][ph] * (*_dfluid_viscosity_dvar)[_qp][ph][pvar] /
     232             :                    std::pow(_fluid_viscosity[_qp][ph], 2)) *
     233             :               _phi[_j][_qp];
     234             : 
     235    20344288 :       Real dvelocity_abs = 0.0;
     236    20344288 :       if (velocity_abs > 0.0)
     237    18202624 :         dvelocity_abs = velocity * dvelocity / velocity_abs;
     238             : 
     239    20344288 :       Real ddiffusion = _phi[_j][_qp] * (*_dporosity_qp_dvar)[_qp][pvar] * _tortuosity[_qp][ph] *
     240    20344288 :                         _diffusion_coeff[_qp][ph][_fluid_component];
     241    20344288 :       ddiffusion += _phi[_j][_qp] * _porosity_qp[_qp] * (*_dtortuosity_dvar)[_qp][ph][pvar] *
     242             :                     _diffusion_coeff[_qp][ph][_fluid_component];
     243    20344288 :       ddiffusion += _phi[_j][_qp] * _porosity_qp[_qp] * _tortuosity[_qp][ph] *
     244    20344288 :                     (*_ddiffusion_coeff_dvar)[_qp][ph][_fluid_component][pvar];
     245    20344288 :       ddiffusion += _disp_trans[ph] * dvelocity_abs;
     246             : 
     247    20344288 :       RankTwoTensor ddispersion;
     248             :       ddispersion.zero();
     249    20344288 :       if (velocity_abs > 0.0)
     250             :       {
     251    18202624 :         RankTwoTensor dv2a, dv2b;
     252    18202624 :         dv2a = RankTwoTensor::outerProduct(velocity, dvelocity);
     253    18202624 :         dv2b = RankTwoTensor::outerProduct(dvelocity, velocity);
     254    18202624 :         ddispersion = (_disp_long[ph] - _disp_trans[ph]) * (dv2a + dv2b) / velocity_abs;
     255    36405248 :         ddispersion -=
     256    18202624 :             (_disp_long[ph] - _disp_trans[ph]) * v2 * dvelocity_abs / velocity_abs / velocity_abs;
     257             :       }
     258             : 
     259    20344288 :       dflux += _phi[_j][_qp] * (*_dfluid_density_qp_dvar)[_qp][ph][pvar] *
     260    20344288 :                (diffusion * _identity_tensor + dispersion) *
     261    20344288 :                _grad_mass_frac[_qp][ph][_fluid_component];
     262    20344288 :       dflux += _fluid_density_qp[_qp][ph] * (ddiffusion * _identity_tensor + ddispersion) *
     263             :                _grad_mass_frac[_qp][ph][_fluid_component];
     264             : 
     265             :       // NOTE: Here we assume that d(grad_mass_frac)/d(var) = d(mass_frac)/d(var) * grad_phi
     266             :       //       This is true for most PorousFlow scenarios, but not for chemical reactions
     267             :       //       where mass_frac is a nonlinear function of the primary MOOSE Variables
     268             :       dflux += _fluid_density_qp[_qp][ph] * (diffusion * _identity_tensor + dispersion) *
     269    20344288 :                (*_dmass_frac_dvar)[_qp][ph][_fluid_component][pvar] * _grad_phi[_j][_qp];
     270             :     }
     271             : 
     272    19229104 :     return _grad_test[_i][_qp] * dflux;
     273             :   }
     274             :   else
     275             :     libmesh_ignore(jvar);
     276           0 :   return 0.0;
     277             : }
     278             : 
     279             : template class PorousFlowDispersiveFluxTempl<false>;
     280             : template class PorousFlowDispersiveFluxTempl<true>;