LCOV - code coverage report
Current view: top level - src/bcs - PorousFlowSink.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 172 185 93.0 %
Date: 2017-11-18 13:30:36 Functions: 10 10 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 "PorousFlowSink.h"
       9             : 
      10             : // MOOSE includes
      11             : #include "MooseVariable.h"
      12             : 
      13             : #include "libmesh/quadrature.h"
      14             : 
      15             : // C++ includes
      16             : #include <iostream>
      17             : 
      18             : template <>
      19             : InputParameters
      20          55 : validParams<PorousFlowSink>()
      21             : {
      22          55 :   InputParameters params = validParams<IntegratedBC>();
      23         165 :   params.addRequiredParam<UserObjectName>(
      24          55 :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      25         165 :   params.addParam<unsigned int>("fluid_phase",
      26             :                                 "If supplied, then this BC will potentially be a function of fluid "
      27             :                                 "pressure, and you can use mass_fraction_component, use_mobility, "
      28             :                                 "use_relperm, use_enthalpy and use_energy.  If not supplied, then "
      29          55 :                                 "this BC can only be a function of temperature");
      30         165 :   params.addParam<unsigned int>("mass_fraction_component",
      31             :                                 "The index corresponding to a fluid "
      32             :                                 "component.  If supplied, the flux will "
      33             :                                 "be multiplied by the nodal mass "
      34          55 :                                 "fraction for the component");
      35         165 :   params.addParam<bool>("use_mobility",
      36             :                         false,
      37             :                         "If true, then fluxes are multiplied by "
      38             :                         "(density*permeability_nn/viscosity), where the "
      39             :                         "'_nn' indicates the component normal to the "
      40             :                         "boundary.  In this case bare_flux is measured in "
      41             :                         "Pa.m^-1.  This can be used in conjunction with "
      42          55 :                         "other use_*");
      43         165 :   params.addParam<bool>("use_relperm",
      44             :                         false,
      45             :                         "If true, then fluxes are multiplied by relative "
      46             :                         "permeability.  This can be used in conjunction with "
      47          55 :                         "other use_*");
      48         165 :   params.addParam<bool>("use_enthalpy",
      49             :                         false,
      50             :                         "If true, then fluxes are multiplied by enthalpy.  "
      51             :                         "In this case bare_flux is measured in kg.m^-2.s^-1 "
      52             :                         "/ (J.kg).  This can be used in conjunction with "
      53          55 :                         "other use_*");
      54         165 :   params.addParam<bool>("use_internal_energy",
      55             :                         false,
      56             :                         "If true, then fluxes are multiplied by fluid internal energy. "
      57             :                         " In this case bare_flux is measured in kg.m^-2.s^-1 / (J.kg). "
      58          55 :                         " This can be used in conjunction with other use_*");
      59         165 :   params.addParam<bool>("use_thermal_conductivity",
      60             :                         false,
      61             :                         "If true, then fluxes are multiplied by "
      62             :                         "thermal conductivity projected onto "
      63             :                         "the normal direction.  This can be "
      64          55 :                         "used in conjunction with other use_*");
      65         165 :   params.addParam<FunctionName>(
      66             :       "flux_function",
      67             :       1.0,
      68             :       "The flux.  The flux is OUT of the medium: hence positive values of "
      69             :       "this function means this BC will act as a SINK, while negative values "
      70             :       "indicate this flux will be a SOURCE.  The functional form is useful "
      71             :       "for spatially or temporally varying sinks.  Without any use_*, this "
      72             :       "function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case "
      73          55 :       "with only heat and no fluids)");
      74         110 :   params.addClassDescription("Applies a flux sink to a boundary.");
      75          55 :   return params;
      76             : }
      77             : 
      78          55 : PorousFlowSink::PorousFlowSink(const InputParameters & parameters)
      79             :   : IntegratedBC(parameters),
      80         110 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      81         110 :     _involves_fluid(isParamValid("fluid_phase")),
      82         206 :     _ph(_involves_fluid ? getParam<unsigned int>("fluid_phase") : 0),
      83         110 :     _use_mass_fraction(isParamValid("mass_fraction_component")),
      84             :     _has_mass_fraction(
      85         159 :         hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
      86         153 :         hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      87             :             "dPorousFlow_mass_frac_nodal_dvar")),
      88         170 :     _sp(_use_mass_fraction ? getParam<unsigned int>("mass_fraction_component") : 0),
      89         165 :     _use_mobility(getParam<bool>("use_mobility")),
      90             :     _has_mobility(
      91         154 :         hasMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp") &&
      92         187 :         hasMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar") &&
      93         187 :         hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
      94         143 :         hasMaterialProperty<std::vector<std::vector<Real>>>(
      95          44 :             "dPorousFlow_fluid_phase_density_nodal_dvar") &&
      96         242 :         hasMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal") &&
      97         143 :         hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
      98         165 :     _use_relperm(getParam<bool>("use_relperm")),
      99         154 :     _has_relperm(hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
     100         143 :                  hasMaterialProperty<std::vector<std::vector<Real>>>(
     101             :                      "dPorousFlow_relative_permeability_nodal_dvar")),
     102         165 :     _use_enthalpy(getParam<bool>("use_enthalpy")),
     103         120 :     _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
     104          75 :                   hasMaterialProperty<std::vector<std::vector<Real>>>(
     105             :                       "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
     106         165 :     _use_internal_energy(getParam<bool>("use_internal_energy")),
     107             :     _has_internal_energy(
     108         120 :         hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_internal_energy_nodal") &&
     109          75 :         hasMaterialProperty<std::vector<std::vector<Real>>>(
     110             :             "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")),
     111         165 :     _use_thermal_conductivity(getParam<bool>("use_thermal_conductivity")),
     112             :     _has_thermal_conductivity(
     113         117 :         hasMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp") &&
     114          69 :         hasMaterialProperty<std::vector<RealTensorValue>>(
     115             :             "dPorousFlow_thermal_conductivity_qp_dvar")),
     116         110 :     _m_func(getFunction("flux_function")),
     117          55 :     _permeability(_has_mobility
     118         154 :                       ? &getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")
     119             :                       : nullptr),
     120             :     _dpermeability_dvar(
     121          55 :         _has_mobility
     122         154 :             ? &getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")
     123             :             : nullptr),
     124          55 :     _dpermeability_dgradvar(_has_mobility
     125          99 :                                 ? &getMaterialProperty<std::vector<std::vector<RealTensorValue>>>(
     126         143 :                                       "dPorousFlow_permeability_qp_dgradvar")
     127             :                                 : nullptr),
     128             :     _fluid_density_node(
     129          55 :         _has_mobility
     130         154 :             ? &getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")
     131             :             : nullptr),
     132          55 :     _dfluid_density_node_dvar(_has_mobility
     133          99 :                                   ? &getMaterialProperty<std::vector<std::vector<Real>>>(
     134         143 :                                         "dPorousFlow_fluid_phase_density_nodal_dvar")
     135             :                                   : nullptr),
     136          55 :     _fluid_viscosity(_has_mobility
     137         154 :                          ? &getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
     138             :                          : nullptr),
     139          55 :     _dfluid_viscosity_dvar(_has_mobility
     140          99 :                                ? &getMaterialProperty<std::vector<std::vector<Real>>>(
     141         143 :                                      "dPorousFlow_viscosity_nodal_dvar")
     142             :                                : nullptr),
     143             :     _relative_permeability(
     144          55 :         _has_relperm
     145         154 :             ? &getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal")
     146             :             : nullptr),
     147          55 :     _drelative_permeability_dvar(_has_relperm
     148          99 :                                      ? &getMaterialProperty<std::vector<std::vector<Real>>>(
     149         143 :                                            "dPorousFlow_relative_permeability_nodal_dvar")
     150             :                                      : nullptr),
     151             :     _mass_fractions(
     152          55 :         _has_mass_fraction
     153         159 :             ? &getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")
     154             :             : nullptr),
     155          55 :     _dmass_fractions_dvar(_has_mass_fraction
     156         104 :                               ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
     157         153 :                                     "dPorousFlow_mass_frac_nodal_dvar")
     158             :                               : nullptr),
     159             :     _enthalpy(
     160          55 :         _has_enthalpy
     161         130 :             ? &getMaterialPropertyByName<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal")
     162             :             : nullptr),
     163          55 :     _denthalpy_dvar(_has_enthalpy
     164          65 :                         ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
     165          65 :                               "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
     166             :                         : nullptr),
     167          55 :     _internal_energy(_has_internal_energy
     168          65 :                          ? &getMaterialPropertyByName<std::vector<Real>>(
     169          65 :                                "PorousFlow_fluid_phase_internal_energy_nodal")
     170             :                          : nullptr),
     171          55 :     _dinternal_energy_dvar(_has_internal_energy
     172          65 :                                ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
     173          65 :                                      "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
     174             :                                : nullptr),
     175             :     _thermal_conductivity(
     176          55 :         _has_thermal_conductivity
     177         117 :             ? &getMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp")
     178             :             : nullptr),
     179          55 :     _dthermal_conductivity_dvar(_has_thermal_conductivity
     180          62 :                                     ? &getMaterialProperty<std::vector<RealTensorValue>>(
     181          69 :                                           "dPorousFlow_thermal_conductivity_qp_dvar")
     182        2105 :                                     : nullptr)
     183             : {
     184          55 :   if (_involves_fluid && _ph >= _dictator.numPhases())
     185           0 :     mooseError("PorousFlowSink: The Dictator declares that the number of fluid phases is ",
     186           0 :                _dictator.numPhases(),
     187             :                ", but you have set the fluid_phase to ",
     188             :                _ph,
     189           0 :                ".  You must try harder.");
     190          62 :   if (!_involves_fluid && (_use_mass_fraction || _use_mobility || _use_relperm || _use_enthalpy ||
     191           7 :                            _use_internal_energy))
     192             :     mooseError("PorousFlowSink: To use_mass_fraction, use_mobility, use_relperm, use_enthalpy or "
     193           0 :                "use_internal_energy, you must provide a fluid phase number");
     194          55 :   if (_use_mass_fraction && _sp >= _dictator.numComponents())
     195           0 :     mooseError("PorousFlowSink: The Dictator declares that the number of fluid components is ",
     196           0 :                _dictator.numComponents(),
     197             :                ", but you have set the mass_fraction_component to ",
     198             :                _sp,
     199           0 :                ".  Please be assured that the Dictator has noted your error.");
     200          55 :   if (_use_mass_fraction && !_has_mass_fraction)
     201             :     mooseError("PorousFlowSink: You have used the use_mass_fraction flag, but you have no "
     202           0 :                "mass_fraction Material");
     203          55 :   if (_use_mobility && !_has_mobility)
     204             :     mooseError("PorousFlowSink: You have used the use_mobility flag, but there are not the "
     205           0 :                "required Materials for this");
     206          55 :   if (_use_relperm && !_has_relperm)
     207             :     mooseError(
     208           0 :         "PorousFlowSink: You have used the use_relperm flag, but you have no relperm Material");
     209          55 :   if (_use_enthalpy && !_has_enthalpy)
     210             :     mooseError(
     211           0 :         "PorousFlowSink: You have used the use_enthalpy flag, but you have no enthalpy Material");
     212          55 :   if (_use_internal_energy && !_has_internal_energy)
     213             :     mooseError("PorousFlowSink: You have used the use_internal_energy flag, but you have no "
     214           0 :                "internal_energy Material");
     215          55 :   if (_use_thermal_conductivity && !_has_thermal_conductivity)
     216             :     mooseError("PorousFlowSink: You have used the use_thermal_conductivity flag, but you have no "
     217           0 :                "thermal_conductivity Material");
     218          55 : }
     219             : 
     220             : Real
     221      544296 : PorousFlowSink::computeQpResidual()
     222             : {
     223     1088592 :   Real flux = _test[_i][_qp] * multiplier();
     224      544296 :   if (_use_mobility)
     225             :   {
     226             :     const Real k =
     227     1145736 :         ((*_permeability)[_qp] * _normals[_qp]) * _normals[_qp]; // do not upwind permeability
     228     1145736 :     flux *= (*_fluid_density_node)[_i][_ph] * k / (*_fluid_viscosity)[_i][_ph];
     229             :   }
     230      544296 :   if (_use_relperm)
     231      778736 :     flux *= (*_relative_permeability)[_i][_ph];
     232      544296 :   if (_use_mass_fraction)
     233      765968 :     flux *= (*_mass_fractions)[_i][_ph][_sp];
     234      544296 :   if (_use_enthalpy)
     235      165856 :     flux *= (*_enthalpy)[_i][_ph];
     236      544296 :   if (_use_internal_energy)
     237       49840 :     flux *= (*_internal_energy)[_i][_ph];
     238      544296 :   if (_use_thermal_conductivity)
     239       49728 :     flux *= ((*_thermal_conductivity)[_qp] * _normals[_qp]) *
     240       24864 :             _normals[_qp]; // do not upwind thermal_conductivity
     241             : 
     242      544296 :   return flux;
     243             : }
     244             : 
     245             : Real
     246      841776 : PorousFlowSink::computeQpJacobian()
     247             : {
     248      841776 :   return jac(_var.number());
     249             : }
     250             : 
     251             : Real
     252      228912 : PorousFlowSink::computeQpOffDiagJacobian(unsigned int jvar)
     253             : {
     254      228912 :   return jac(jvar);
     255             : }
     256             : 
     257             : Real
     258     1070688 : PorousFlowSink::jac(unsigned int jvar) const
     259             : {
     260     1070688 :   if (_dictator.notPorousFlowVariable(jvar))
     261             :     return 0.0;
     262     1070688 :   const unsigned int pvar = _dictator.porousFlowVariableNum(jvar);
     263             : 
     264             :   // For _i != _j, note:
     265             :   // since the only non-upwinded contribution to the residual is
     266             :   // from the permeability and thermal_conductivity, the only contribution
     267             :   // of the residual at node _i from changing jvar at node _j is through
     268             :   // the derivative of permeability or thermal_conductivity
     269             : 
     270     2141376 :   Real flux = _test[_i][_qp] * multiplier();
     271     2141376 :   Real deriv = _test[_i][_qp] * (_i != _j ? 0.0 : dmultiplier_dvar(pvar));
     272             : 
     273     1070688 :   if (_use_mobility)
     274             :   {
     275      764448 :     const Real k = ((*_permeability)[_qp] * _normals[_qp]) * _normals[_qp];
     276      764448 :     const Real mob = (*_fluid_density_node)[_i][_ph] * k / (*_fluid_viscosity)[_i][_ph];
     277      509632 :     RealTensorValue ktprime = (*_dpermeability_dvar)[_qp][pvar] * _phi[_j][_qp];
     278     1783712 :     for (unsigned i = 0; i < LIBMESH_DIM; ++i)
     279     2293344 :       ktprime += (*_dpermeability_dgradvar)[_qp][i][pvar] * _grad_phi[_j][_qp](i);
     280      254816 :     const Real kprime = (ktprime * _normals[_qp]) * _normals[_qp];
     281             : 
     282      254816 :     Real mobprime = (*_fluid_density_node)[_i][_ph] * kprime / (*_fluid_viscosity)[_i][_ph];
     283      254816 :     mobprime +=
     284             :         (_i != _j
     285      286744 :              ? 0.0
     286       63856 :              : (*_dfluid_density_node_dvar)[_i][_ph][pvar] * k / (*_fluid_viscosity)[_i][_ph] -
     287       95784 :                    (*_fluid_density_node)[_i][_ph] * k * (*_dfluid_viscosity_dvar)[_i][_ph][pvar] /
     288             :                        std::pow((*_fluid_viscosity)[_i][_ph], 2));
     289      254816 :     deriv = mob * deriv + mobprime * flux;
     290      254816 :     flux *= mob;
     291             :   }
     292     1070688 :   if (_use_relperm)
     293             :   {
     294      355576 :     const Real relperm_prime = (_i != _j ? 0.0 : (*_drelative_permeability_dvar)[_i][_ph][pvar]);
     295      632000 :     deriv = (*_relative_permeability)[_i][_ph] * deriv + relperm_prime * flux;
     296      316000 :     flux *= (*_relative_permeability)[_i][_ph];
     297             :   }
     298     1070688 :   if (_use_mass_fraction)
     299             :   {
     300      330576 :     const Real mf_prime = (_i != _j ? 0.0 : (*_dmass_fractions_dvar)[_i][_ph][_sp][pvar]);
     301      585536 :     deriv = (*_mass_fractions)[_i][_ph][_sp] * deriv + mf_prime * flux;
     302      292768 :     flux *= (*_mass_fractions)[_i][_ph][_sp];
     303             :   }
     304     1070688 :   if (_use_enthalpy)
     305             :   {
     306       60800 :     const Real en_prime = (_i != _j ? 0.0 : (*_denthalpy_dvar)[_i][_ph][pvar]);
     307      108032 :     deriv = (*_enthalpy)[_i][_ph] * deriv + en_prime * flux;
     308       54016 :     flux *= (*_enthalpy)[_i][_ph];
     309             :   }
     310     1070688 :   if (_use_internal_energy)
     311             :   {
     312       18544 :     const Real ie_prime = (_i != _j ? 0.0 : (*_dinternal_energy_dvar)[_i][_ph][pvar]);
     313       32896 :     deriv = (*_internal_energy)[_i][_ph] * deriv + ie_prime * flux;
     314       16448 :     flux *= (*_internal_energy)[_i][_ph];
     315             :   }
     316     1070688 :   if (_use_thermal_conductivity)
     317             :   {
     318       48384 :     const Real tc = ((*_thermal_conductivity)[_qp] * _normals[_qp]) * _normals[_qp];
     319       32256 :     const RealTensorValue tctprime = (*_dthermal_conductivity_dvar)[_qp][pvar] * _phi[_j][_qp];
     320       16128 :     const Real tcprime = (tctprime * _normals[_qp]) * _normals[_qp];
     321       16128 :     deriv = tc * deriv + tcprime * flux;
     322             :     // don't need this: flux *= tc;
     323             :   }
     324             :   return deriv;
     325             : }
     326             : 
     327             : Real
     328     1569132 : PorousFlowSink::multiplier() const
     329             : {
     330     3138264 :   return _m_func.value(_t, _q_point[_qp]);
     331             : }
     332             : 
     333             : Real
     334      130452 : PorousFlowSink::dmultiplier_dvar(unsigned int /*pvar*/) const
     335             : {
     336      130452 :   return 0.0;
     337        2499 : }

Generated by: LCOV version 1.11