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 "PorousFlowDarcyVelocityMaterial.h"
      11             : 
      12             : registerMooseObject("PorousFlowApp", PorousFlowDarcyVelocityMaterial);
      13             : registerMooseObject("PorousFlowApp", ADPorousFlowDarcyVelocityMaterial);
      14             : 
      15             : template <bool is_ad>
      16             : InputParameters
      17         280 : PorousFlowDarcyVelocityMaterialTempl<is_ad>::validParams()
      18             : {
      19         280 :   InputParameters params = PorousFlowMaterial::validParams();
      20         560 :   params.addRequiredParam<RealVectorValue>("gravity",
      21             :                                            "Gravitational acceleration vector downwards (m/s^2)");
      22         280 :   params.addClassDescription("This Material calculates the Darcy velocity for all phases");
      23         560 :   params.addPrivateParam<std::string>("pf_material_type", "darcy_velocity");
      24         280 :   params.set<bool>("at_nodes") = false;
      25         280 :   return params;
      26           0 : }
      27             : 
      28             : template <bool is_ad>
      29         218 : PorousFlowDarcyVelocityMaterialTempl<is_ad>::PorousFlowDarcyVelocityMaterialTempl(
      30             :     const InputParameters & parameters)
      31             :   : PorousFlowMaterial(parameters),
      32         436 :     _num_phases(_dictator.numPhases()),
      33         218 :     _num_var(_dictator.numVariables()),
      34         218 :     _permeability(getGenericMaterialProperty<RealTensorValue, is_ad>("PorousFlow_permeability_qp")),
      35         218 :     _dpermeability_dvar(is_ad ? nullptr
      36         197 :                               : &getMaterialProperty<std::vector<RealTensorValue>>(
      37             :                                     "dPorousFlow_permeability_qp_dvar")),
      38         218 :     _dpermeability_dgradvar(is_ad ? nullptr
      39         197 :                                   : &getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
      40             :                                         "dPorousFlow_permeability_qp_dgradvar")),
      41         218 :     _fluid_density(
      42         218 :         getGenericMaterialProperty<std::vector<Real>, is_ad>("PorousFlow_fluid_phase_density_qp")),
      43         218 :     _dfluid_density_dvar(is_ad ? nullptr
      44         197 :                                : &getMaterialProperty<std::vector<std::vector<Real>>>(
      45             :                                      "dPorousFlow_fluid_phase_density_qp_dvar")),
      46         218 :     _fluid_viscosity(
      47         218 :         getGenericMaterialProperty<std::vector<Real>, is_ad>("PorousFlow_viscosity_qp")),
      48         218 :     _dfluid_viscosity_dvar(is_ad ? nullptr
      49         197 :                                  : &getMaterialProperty<std::vector<std::vector<Real>>>(
      50             :                                        "dPorousFlow_viscosity_qp_dvar")),
      51         218 :     _relative_permeability(getGenericMaterialProperty<std::vector<Real>, is_ad>(
      52             :         "PorousFlow_relative_permeability_qp")),
      53         218 :     _drelative_permeability_dvar(is_ad ? nullptr
      54         197 :                                        : &getMaterialProperty<std::vector<std::vector<Real>>>(
      55             :                                              "dPorousFlow_relative_permeability_qp_dvar")),
      56         218 :     _grad_p(getGenericMaterialProperty<std::vector<RealGradient>, is_ad>(
      57             :         "PorousFlow_grad_porepressure_qp")),
      58         218 :     _dgrad_p_dgradvar(is_ad ? nullptr
      59         197 :                             : &getMaterialProperty<std::vector<std::vector<Real>>>(
      60             :                                   "dPorousFlow_grad_porepressure_qp_dgradvar")),
      61         218 :     _dgrad_p_dvar(is_ad ? nullptr
      62         197 :                         : &getMaterialProperty<std::vector<std::vector<RealGradient>>>(
      63             :                               "dPorousFlow_grad_porepressure_qp_dvar")),
      64         436 :     _gravity(getParam<RealVectorValue>("gravity")),
      65         436 :     _darcy_velocity(declareGenericProperty<std::vector<RealVectorValue>, is_ad>(
      66             :         "PorousFlow_darcy_velocity_qp")),
      67         218 :     _ddarcy_velocity_dvar(is_ad ? nullptr
      68         197 :                                 : &declareProperty<std::vector<std::vector<RealVectorValue>>>(
      69             :                                       "dPorousFlow_darcy_velocity_qp_dvar")),
      70         218 :     _ddarcy_velocity_dgradvar(
      71             :         is_ad ? nullptr
      72         197 :               : &declareProperty<std::vector<std::vector<std::vector<RealVectorValue>>>>(
      73         218 :                     "dPorousFlow_darcy_velocity_qp_dgradvar"))
      74             : {
      75         218 :   if (_nodal_material == true)
      76           2 :     mooseError("PorousFlowDarcyVelocityMaterial is only defined for at_nodes = false");
      77         216 : }
      78             : 
      79             : template <bool is_ad>
      80             : void
      81        4240 : PorousFlowDarcyVelocityMaterialTempl<is_ad>::computeQpProperties()
      82             : {
      83        4240 :   _darcy_velocity[_qp].resize(_num_phases);
      84             : 
      85       10428 :   for (unsigned ph = 0; ph < _num_phases; ++ph)
      86        6260 :     _darcy_velocity[_qp][ph] =
      87        6188 :         -(_permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density[_qp][ph] * _gravity) *
      88        6188 :           _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph]);
      89             : 
      90             :   if constexpr (!is_ad)
      91             :   {
      92        4168 :     (*_ddarcy_velocity_dvar)[_qp].resize(_num_phases);
      93       10284 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
      94        6116 :       (*_ddarcy_velocity_dvar)[_qp][ph].resize(_num_var);
      95             : 
      96       10284 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
      97        6860 :       for (unsigned v = 0; v < _num_var; ++v)
      98             :       {
      99         744 :         (*_ddarcy_velocity_dvar)[_qp][ph][v] =
     100         744 :             -(*_dpermeability_dvar)[_qp][v] *
     101         744 :             (_grad_p[_qp][ph] - _fluid_density[_qp][ph] * _gravity) *
     102         744 :             _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph];
     103             :         (*_ddarcy_velocity_dvar)[_qp][ph][v] -=
     104        1488 :             _permeability[_qp] *
     105         744 :             ((*_dgrad_p_dvar)[_qp][ph][v] - (*_dfluid_density_dvar)[_qp][ph][v] * _gravity) *
     106         744 :             _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph];
     107             :         (*_ddarcy_velocity_dvar)[_qp][ph][v] -=
     108        1488 :             _permeability[_qp] * (_grad_p[_qp][ph] - _fluid_density[_qp][ph] * _gravity) *
     109         744 :             ((*_drelative_permeability_dvar)[_qp][ph][v] / _fluid_viscosity[_qp][ph] -
     110         744 :              _relative_permeability[_qp][ph] * (*_dfluid_viscosity_dvar)[_qp][ph][v] /
     111             :                  std::pow(_fluid_viscosity[_qp][ph], 2));
     112             :       }
     113             : 
     114        4168 :     (*_ddarcy_velocity_dgradvar)[_qp].resize(_num_phases);
     115       10284 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     116             :     {
     117        6116 :       (*_ddarcy_velocity_dgradvar)[_qp][ph].resize(LIBMESH_DIM);
     118       24464 :       for (unsigned i = 0; i < LIBMESH_DIM; ++i)
     119       18348 :         (*_ddarcy_velocity_dgradvar)[_qp][ph][i].resize(_num_var);
     120             :     }
     121             : 
     122       10284 :     for (unsigned ph = 0; ph < _num_phases; ++ph)
     123       24464 :       for (unsigned i = 0; i < LIBMESH_DIM; ++i)
     124       20580 :         for (unsigned v = 0; v < _num_var; ++v)
     125             :         {
     126        2232 :           (*_ddarcy_velocity_dgradvar)[_qp][ph][i][v] =
     127        2232 :               -(*_dpermeability_dgradvar)[_qp][i][v] *
     128        2232 :               (_grad_p[_qp][ph] - _fluid_density[_qp][ph] * _gravity);
     129        8928 :           for (unsigned j = 0; j < LIBMESH_DIM; ++j)
     130        6696 :             (*_ddarcy_velocity_dgradvar)[_qp][ph][i][v](j) -=
     131        6696 :                 _permeability[_qp](j, i) * (*_dgrad_p_dgradvar)[_qp][ph][v];
     132             :           (*_ddarcy_velocity_dgradvar)[_qp][ph][i][v] *=
     133        2232 :               _relative_permeability[_qp][ph] / _fluid_viscosity[_qp][ph];
     134             :         }
     135             :   }
     136        4240 : }
     137             : 
     138             : template class PorousFlowDarcyVelocityMaterialTempl<false>;
     139             : template class PorousFlowDarcyVelocityMaterialTempl<true>;