LCOV - code coverage report
Current view: top level - src/dirackernels - PorousFlowLineSink.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 190 190 100.0 %
Date: 2017-11-18 13:30:36 Functions: 11 11 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 "PorousFlowLineSink.h"
       9             : #include "libmesh/utility.h"
      10             : 
      11             : template <>
      12             : InputParameters
      13          62 : validParams<PorousFlowLineSink>()
      14             : {
      15          62 :   InputParameters params = validParams<PorousFlowLineGeometry>();
      16         186 :   MooseEnum p_or_t_choice("pressure=0 temperature=1", "pressure");
      17         186 :   params.addParam<MooseEnum>("function_of",
      18             :                              p_or_t_choice,
      19             :                              "Modifying functions will be a function of either pressure and "
      20             :                              "permeability (eg, for boreholes that pump fluids) or "
      21             :                              "temperature and thermal conductivity (eg, for boreholes that "
      22          62 :                              "pump pure heat with no fluid flow)");
      23         186 :   params.addRequiredParam<UserObjectName>(
      24             :       "SumQuantityUO",
      25             :       "User Object of type=PorousFlowSumQuantity in which to place the total "
      26          62 :       "outflow from the line sink for each time step.");
      27         186 :   params.addRequiredParam<UserObjectName>(
      28          62 :       "PorousFlowDictator", "The UserObject that holds the list of PorousFlow variable names");
      29         186 :   params.addParam<unsigned int>(
      30             :       "fluid_phase",
      31             :       0,
      32             :       "The fluid phase whose pressure (and potentially mobility, enthalpy, etc) "
      33             :       "controls the flux to the line sink.  For p_or_t=temperature, and without "
      34          62 :       "any use_*, this parameter is irrelevant");
      35         186 :   params.addParam<unsigned int>("mass_fraction_component",
      36             :                                 "The index corresponding to a fluid "
      37             :                                 "component.  If supplied, the flux will "
      38             :                                 "be multiplied by the nodal mass "
      39          62 :                                 "fraction for the component");
      40         186 :   params.addParam<bool>(
      41          62 :       "use_relative_permeability", false, "Multiply the flux by the fluid relative permeability");
      42         186 :   params.addParam<bool>("use_mobility", false, "Multiply the flux by the fluid mobility");
      43         186 :   params.addParam<bool>("use_enthalpy", false, "Multiply the flux by the fluid enthalpy");
      44         186 :   params.addParam<bool>(
      45          62 :       "use_internal_energy", false, "Multiply the flux by the fluid internal energy");
      46         124 :   params.addClassDescription("Approximates a line sink in the mesh by a sequence of weighted Dirac "
      47          62 :                              "points whose positions are read from a file");
      48          62 :   return params;
      49             : }
      50             : 
      51          62 : PorousFlowLineSink::PorousFlowLineSink(const InputParameters & parameters)
      52             :   : PorousFlowLineGeometry(parameters),
      53         120 :     _dictator(getUserObject<PorousFlowDictator>("PorousFlowDictator")),
      54             :     _total_outflow_mass(
      55         120 :         const_cast<PorousFlowSumQuantity &>(getUserObject<PorousFlowSumQuantity>("SumQuantityUO"))),
      56             : 
      57             :     _has_porepressure(
      58         178 :         hasMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp") &&
      59         176 :         hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_porepressure_qp_dvar")),
      60         154 :     _has_temperature(hasMaterialProperty<Real>("PorousFlow_temperature_qp") &&
      61         128 :                      hasMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
      62             :     _has_mass_fraction(
      63         180 :         hasMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal") &&
      64         180 :         hasMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
      65             :             "dPorousFlow_mass_frac_nodal_dvar")),
      66             :     _has_relative_permeability(
      67         174 :         hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
      68         168 :         hasMaterialProperty<std::vector<std::vector<Real>>>(
      69             :             "dPorousFlow_relative_permeability_nodal_dvar")),
      70             :     _has_mobility(
      71         174 :         hasMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal") &&
      72         168 :         hasMaterialProperty<std::vector<std::vector<Real>>>(
      73          54 :             "dPorousFlow_relative_permeability_nodal_dvar") &&
      74         222 :         hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal") &&
      75         168 :         hasMaterialProperty<std::vector<std::vector<Real>>>(
      76          53 :             "dPorousFlow_fluid_phase_density_nodal_dvar") &&
      77         274 :         hasMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal") &&
      78         156 :         hasMaterialProperty<std::vector<std::vector<Real>>>("dPorousFlow_viscosity_nodal_dvar")),
      79         150 :     _has_enthalpy(hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal") &&
      80         120 :                   hasMaterialProperty<std::vector<std::vector<Real>>>(
      81             :                       "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")),
      82             :     _has_internal_energy(
      83         150 :         hasMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_internal_energy_nodal") &&
      84         120 :         hasMaterialProperty<std::vector<std::vector<Real>>>(
      85             :             "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")),
      86             : 
      87         120 :     _p_or_t(getParam<MooseEnum>("function_of").getEnum<PorTchoice>()),
      88         120 :     _use_mass_fraction(isParamValid("mass_fraction_component")),
      89         180 :     _use_relative_permeability(getParam<bool>("use_relative_permeability")),
      90         180 :     _use_mobility(getParam<bool>("use_mobility")),
      91         180 :     _use_enthalpy(getParam<bool>("use_enthalpy")),
      92         180 :     _use_internal_energy(getParam<bool>("use_internal_energy")),
      93             : 
      94         180 :     _ph(getParam<unsigned int>("fluid_phase")),
      95         156 :     _sp(_use_mass_fraction ? getParam<unsigned int>("mass_fraction_component") : 0),
      96             : 
      97         110 :     _pp((_p_or_t == PorTchoice::pressure && _has_porepressure)
      98         218 :             ? &getMaterialProperty<std::vector<Real>>("PorousFlow_porepressure_qp")
      99             :             : nullptr),
     100         110 :     _dpp_dvar((_p_or_t == PorTchoice::pressure && _has_porepressure)
     101         109 :                   ? &getMaterialProperty<std::vector<std::vector<Real>>>(
     102         158 :                         "dPorousFlow_porepressure_qp_dvar")
     103             :                   : nullptr),
     104          70 :     _temperature((_p_or_t == PorTchoice::temperature && _has_temperature)
     105         138 :                      ? &getMaterialProperty<Real>("PorousFlow_temperature_qp")
     106             :                      : nullptr),
     107             :     _dtemperature_dvar(
     108          70 :         (_p_or_t == PorTchoice::temperature && _has_temperature)
     109         138 :             ? &getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")
     110             :             : nullptr),
     111             :     _fluid_density_node(
     112          80 :         (_use_mobility && _has_mobility)
     113         154 :             ? &getMaterialProperty<std::vector<Real>>("PorousFlow_fluid_phase_density_nodal")
     114             :             : nullptr),
     115          80 :     _dfluid_density_node_dvar((_use_mobility && _has_mobility)
     116          77 :                                   ? &getMaterialProperty<std::vector<std::vector<Real>>>(
     117          94 :                                         "dPorousFlow_fluid_phase_density_nodal_dvar")
     118             :                                   : nullptr),
     119          80 :     _fluid_viscosity((_use_mobility && _has_mobility)
     120         154 :                          ? &getMaterialProperty<std::vector<Real>>("PorousFlow_viscosity_nodal")
     121             :                          : nullptr),
     122          80 :     _dfluid_viscosity_dvar((_use_mobility && _has_mobility)
     123          77 :                                ? &getMaterialProperty<std::vector<std::vector<Real>>>(
     124          94 :                                      "dPorousFlow_viscosity_nodal_dvar")
     125             :                                : nullptr),
     126             :     _relative_permeability(
     127         123 :         ((_use_mobility && _has_mobility) ||
     128          53 :          (_use_relative_permeability && _has_relative_permeability))
     129         172 :             ? &getMaterialProperty<std::vector<Real>>("PorousFlow_relative_permeability_nodal")
     130             :             : nullptr),
     131         123 :     _drelative_permeability_dvar(((_use_mobility && _has_mobility) ||
     132          53 :                                   (_use_relative_permeability && _has_relative_permeability))
     133          86 :                                      ? &getMaterialProperty<std::vector<std::vector<Real>>>(
     134         112 :                                            "dPorousFlow_relative_permeability_nodal_dvar")
     135             :                                      : nullptr),
     136             :     _mass_fractions(
     137          78 :         (_use_mass_fraction && _has_mass_fraction)
     138         154 :             ? &getMaterialProperty<std::vector<std::vector<Real>>>("PorousFlow_mass_frac_nodal")
     139             :             : nullptr),
     140          78 :     _dmass_fractions_dvar((_use_mass_fraction && _has_mass_fraction)
     141          77 :                               ? &getMaterialProperty<std::vector<std::vector<std::vector<Real>>>>(
     142          94 :                                     "dPorousFlow_mass_frac_nodal_dvar")
     143             :                               : nullptr),
     144             :     _enthalpy(
     145          60 :         _has_enthalpy
     146         180 :             ? &getMaterialPropertyByName<std::vector<Real>>("PorousFlow_fluid_phase_enthalpy_nodal")
     147             :             : nullptr),
     148          60 :     _denthalpy_dvar(_has_enthalpy
     149          90 :                         ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
     150          90 :                               "dPorousFlow_fluid_phase_enthalpy_nodal_dvar")
     151             :                         : nullptr),
     152          60 :     _internal_energy(_has_internal_energy
     153          90 :                          ? &getMaterialPropertyByName<std::vector<Real>>(
     154          90 :                                "PorousFlow_fluid_phase_internal_energy_nodal")
     155             :                          : nullptr),
     156          60 :     _dinternal_energy_dvar(_has_internal_energy
     157          90 :                                ? &getMaterialPropertyByName<std::vector<std::vector<Real>>>(
     158          90 :                                      "dPorousFlow_fluid_phase_internal_energy_nodal_dvar")
     159        2162 :                                : nullptr)
     160             : {
     161             :   // zero the outflow mass
     162          60 :   _total_outflow_mass.zero();
     163             : 
     164          60 :   if (_ph >= _dictator.numPhases())
     165           1 :     mooseError("PorousFlowLineSink: The Dictator declares that the number of fluid phases is ",
     166           1 :                _dictator.numPhases(),
     167             :                ", but you have set the fluid_phase to ",
     168             :                _ph,
     169           1 :                ".  You must try harder.");
     170          59 :   if (_use_mass_fraction && _sp >= _dictator.numComponents())
     171           1 :     mooseError("PorousFlowLineSink: The Dictator declares that the number of fluid components is ",
     172           1 :                _dictator.numComponents(),
     173             :                ", but you have set the mass_fraction_component to ",
     174             :                _sp,
     175           1 :                ".  Please be assured that the Dictator has noted your error.");
     176          58 :   if (_p_or_t == PorTchoice::pressure && !_has_porepressure)
     177             :     mooseError("PorousFlowLineSink: You have specified function_of=porepressure, but you do not "
     178           1 :                "have a quadpoint porepressure material");
     179          57 :   if (_p_or_t == PorTchoice::temperature && !_has_temperature)
     180             :     mooseError("PorousFlowLineSink: You have specified function_of=temperature, but you do not "
     181           1 :                "have a quadpoint temperature material");
     182          56 :   if (_use_mass_fraction && !_has_mass_fraction)
     183             :     mooseError("PorousFlowLineSink: You have specified a fluid component, but do not have a nodal "
     184           1 :                "mass-fraction material");
     185          55 :   if (_use_relative_permeability && !_has_relative_permeability)
     186             :     mooseError("PorousFlowLineSink: You have set use_relative_permeability=true, but do not have a "
     187           1 :                "nodal relative permeability material");
     188          54 :   if (_use_mobility && !_has_mobility)
     189             :     mooseError("PorousFlowLineSink: You have set use_mobility=true, but do not have nodal density, "
     190           3 :                "relative permeability or viscosity material");
     191          51 :   if (_use_enthalpy && !_has_enthalpy)
     192             :     mooseError("PorousFlowLineSink: You have set use_enthalpy=true, but do not have a nodal "
     193           1 :                "enthalpy material");
     194          50 :   if (_use_internal_energy && !_has_internal_energy)
     195             :     mooseError("PorousFlowLineSink: You have set use_internal_energy=true, but do not have a nodal "
     196           1 :                "internal-energy material");
     197             : 
     198             :   // To correctly compute the Jacobian terms,
     199             :   // tell MOOSE that this DiracKernel depends on all the PorousFlow Variables
     200          49 :   const std::vector<MooseVariable *> & coupled_vars = _dictator.getCoupledMooseVars();
     201         785 :   for (unsigned int i = 0; i < coupled_vars.size(); i++)
     202         229 :     addMooseVariableDependency(coupled_vars[i]);
     203          49 : }
     204             : 
     205             : void
     206        7641 : PorousFlowLineSink::addPoints()
     207             : {
     208             :   // This function gets called just before the DiracKernel is evaluated
     209             :   // so this is a handy place to zero this out.
     210        7641 :   _total_outflow_mass.zero();
     211             : 
     212        7641 :   PorousFlowLineGeometry::addPoints();
     213        7641 : }
     214             : 
     215             : Real
     216      104021 : PorousFlowLineSink::computeQpResidual()
     217             : {
     218             :   // Get the ID we initially assigned to this point
     219      104021 :   const unsigned current_dirac_ptid = currentPointCachedID();
     220      104021 :   Real outflow = computeQpBaseOutflow(current_dirac_ptid);
     221      104020 :   if (outflow == 0.0)
     222             :     return 0.0;
     223             : 
     224       84540 :   if (_use_relative_permeability)
     225       32736 :     outflow *= (*_relative_permeability)[_i][_ph];
     226             : 
     227       84540 :   if (_use_mobility)
     228      107184 :     outflow *= (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] /
     229       26796 :                (*_fluid_viscosity)[_i][_ph];
     230             : 
     231       84540 :   if (_use_mass_fraction)
     232      106120 :     outflow *= (*_mass_fractions)[_i][_ph][_sp];
     233             : 
     234       84540 :   if (_use_enthalpy)
     235       13888 :     outflow *= (*_enthalpy)[_i][_ph];
     236             : 
     237       84540 :   if (_use_internal_energy)
     238       92232 :     outflow *= (*_internal_energy)[_i][_ph];
     239             : 
     240       84540 :   _total_outflow_mass.add(
     241      169080 :       outflow * _dt); // this is not thread safe, but DiracKernel's aren't currently threaded
     242             : 
     243       84540 :   return outflow;
     244             : }
     245             : 
     246             : Real
     247       59312 : PorousFlowLineSink::computeQpJacobian()
     248             : {
     249       59312 :   return jac(_var.number());
     250             : }
     251             : 
     252             : Real
     253       62208 : PorousFlowLineSink::computeQpOffDiagJacobian(unsigned int jvar)
     254             : {
     255       62208 :   return jac(jvar);
     256             : }
     257             : 
     258             : Real
     259      121520 : PorousFlowLineSink::jac(unsigned int jvar)
     260             : {
     261      121520 :   if (_dictator.notPorousFlowVariable(jvar))
     262             :     return 0.0;
     263      121520 :   const unsigned pvar = _dictator.porousFlowVariableNum(jvar);
     264             : 
     265             :   Real outflow;
     266             :   Real outflowp;
     267      121520 :   const unsigned current_dirac_ptid = currentPointCachedID();
     268      121520 :   computeQpBaseOutflowJacobian(jvar, current_dirac_ptid, outflow, outflowp);
     269      121520 :   if (outflow == 0.0 && outflowp == 0.0)
     270             :     return 0.0;
     271             : 
     272      103424 :   if (_use_relative_permeability)
     273             :   {
     274       13152 :     const Real relperm_prime = (_i != _j ? 0.0 : (*_drelative_permeability_dvar)[_i][_ph][pvar]);
     275       23328 :     outflowp = (*_relative_permeability)[_i][_ph] * outflowp + relperm_prime * outflow;
     276       11664 :     outflow *= (*_relative_permeability)[_i][_ph];
     277             :   }
     278             : 
     279      103424 :   if (_use_mobility)
     280             :   {
     281      235136 :     const Real mob = (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] /
     282      117568 :                      (*_fluid_viscosity)[_i][_ph];
     283             :     const Real mob_prime =
     284       58784 :         (_i != _j
     285       66132 :              ? 0.0
     286       14696 :              : (*_drelative_permeability_dvar)[_i][_ph][pvar] * (*_fluid_density_node)[_i][_ph] /
     287        7348 :                        (*_fluid_viscosity)[_i][_ph] +
     288        7348 :                    (*_relative_permeability)[_i][_ph] *
     289       14696 :                        (*_dfluid_density_node_dvar)[_i][_ph][pvar] / (*_fluid_viscosity)[_i][_ph] -
     290        7348 :                    (*_relative_permeability)[_i][_ph] * (*_fluid_density_node)[_i][_ph] *
     291       14696 :                        (*_dfluid_viscosity_dvar)[_i][_ph][pvar] /
     292             :                        Utility::pow<2>((*_fluid_viscosity)[_i][_ph]));
     293       58784 :     outflowp = mob * outflowp + mob_prime * outflow;
     294       58784 :     outflow *= mob;
     295             :   }
     296             : 
     297      103424 :   if (_use_mass_fraction)
     298             :   {
     299             :     const Real mass_fractions_prime =
     300       44604 :         (_i != _j ? 0.0 : (*_dmass_fractions_dvar)[_i][_ph][_sp][pvar]);
     301       79296 :     outflowp = (*_mass_fractions)[_i][_ph][_sp] * outflowp + mass_fractions_prime * outflow;
     302       39648 :     outflow *= (*_mass_fractions)[_i][_ph][_sp];
     303             :   }
     304             : 
     305      103424 :   if (_use_enthalpy)
     306             :   {
     307        6048 :     const Real enthalpy_prime = (_i != _j ? 0.0 : (*_denthalpy_dvar)[_i][_ph][pvar]);
     308       10752 :     outflowp = (*_enthalpy)[_i][_ph] * outflowp + enthalpy_prime * outflow;
     309        5376 :     outflow *= (*_enthalpy)[_i][_ph];
     310             :   }
     311             : 
     312      103424 :   if (_use_internal_energy)
     313             :   {
     314       38556 :     const Real internal_energy_prime = (_i != _j ? 0.0 : (*_dinternal_energy_dvar)[_i][_ph][pvar]);
     315       68544 :     outflowp = (*_internal_energy)[_i][_ph] * outflowp + internal_energy_prime * outflow;
     316             :     // this multiplication was performed, but the code does not need to know: outflow *=
     317             :     // (*_internal_energy)[_i][_ph];
     318             :   }
     319             : 
     320      103424 :   return outflowp;
     321             : }
     322             : 
     323             : Real
     324      244965 : PorousFlowLineSink::ptqp() const
     325             : {
     326      489930 :   return (_p_or_t == PorTchoice::pressure ? (*_pp)[_qp][_ph] : (*_temperature)[_qp]);
     327             : }
     328             : 
     329             : Real
     330      113456 : PorousFlowLineSink::dptqp(unsigned pvar) const
     331             : {
     332      226912 :   return (_p_or_t == PorTchoice::pressure ? (*_dpp_dvar)[_qp][_ph][pvar]
     333      145712 :                                           : (*_dtemperature_dvar)[_qp][pvar]);
     334        2499 : }

Generated by: LCOV version 1.11