LCOV - code coverage report
Current view: top level - src/dirackernels - PorousFlowPeacemanBorehole.C (source / functions) Hit Total Coverage
Test: porous_flow Test Coverage Lines: 131 136 96.3 %
Date: 2017-11-21 14:47:27 Functions: 7 7 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 "PorousFlowPeacemanBorehole.h"
       9             : #include "RotationMatrix.h"
      10             : #include "Function.h"
      11             : 
      12             : template <>
      13             : InputParameters
      14          38 : validParams<PorousFlowPeacemanBorehole>()
      15             : {
      16          38 :   InputParameters params = validParams<PorousFlowLineSink>();
      17         114 :   params.addRequiredParam<FunctionName>(
      18             :       "character",
      19             :       "If zero then borehole does nothing.  If positive the borehole acts as a sink "
      20             :       "(production well) for porepressure > borehole pressure, and does nothing "
      21             :       "otherwise.  If negative the borehole acts as a source (injection well) for "
      22             :       "porepressure < borehole pressure, and does nothing otherwise.  The flow rate "
      23             :       "to/from the borehole is multiplied by |character|, so usually character = +/- "
      24             :       "1, but you can specify other quantities to provide an overall scaling to the "
      25          38 :       "flow if you like.");
      26         114 :   params.addRequiredParam<Real>("bottom_p_or_t",
      27             :                                 "For function_of=pressure, this parameter is the "
      28             :                                 "pressure at the bottom of the borehole, "
      29             :                                 "otherwise it is the temperature at the bottom of "
      30          38 :                                 "the borehole");
      31         114 :   params.addRequiredParam<RealVectorValue>(
      32             :       "unit_weight",
      33             :       "(fluid_density*gravitational_acceleration) as a vector pointing downwards.  "
      34             :       "Note that the borehole pressure at a given z position is bottom_p_or_t + "
      35             :       "unit_weight*(q - q_bottom), where q=(x,y,z) and q_bottom=(x,y,z) of the "
      36             :       "bottom point of the borehole.  The analogous formula holds for "
      37             :       "function_of=temperature.  If you don't want bottomhole pressure (or "
      38             :       "temperature) to vary in the borehole just set unit_weight=0.  Typical value "
      39          38 :       "is = (0,0,-1E4), for water");
      40         114 :   params.addParam<Real>("re_constant",
      41             :                         0.28,
      42             :                         "The dimensionless constant used in evaluating the borehole effective "
      43             :                         "radius.  This depends on the meshing scheme.  Peacemann "
      44             :                         "finite-difference calculations give 0.28, while for rectangular finite "
      45             :                         "elements the result is closer to 0.1594.  (See  Eqn(4.13) of Z Chen, Y "
      46             :                         "Zhang, Well flow models for various numerical methods, Int J Num "
      47          38 :                         "Analysis and Modeling, 3 (2008) 375-388.)");
      48         114 :   params.addParam<Real>("well_constant",
      49             :                         -1.0,
      50             :                         "Usually this is calculated internally from the element geometry, the "
      51             :                         "local borehole direction and segment length, and the permeability.  "
      52             :                         "However, if this parameter is given as a positive number then this "
      53             :                         "number is used instead of the internal calculation.  This speeds up "
      54          38 :                         "computation marginally.  re_constant becomes irrelevant");
      55          76 :   params.addClassDescription("Approximates a borehole in the mesh using the Peaceman approach, ie "
      56             :                              "using a number of point sinks with given radii whose positions are "
      57          38 :                              "read from a file");
      58          38 :   return params;
      59             : }
      60             : 
      61          38 : PorousFlowPeacemanBorehole::PorousFlowPeacemanBorehole(const InputParameters & parameters)
      62             :   : PorousFlowLineSink(parameters),
      63          50 :     _character(getFunction("character")),
      64          75 :     _p_bot(getParam<Real>("bottom_p_or_t")),
      65          50 :     _unit_weight(getParam<RealVectorValue>("unit_weight")),
      66          75 :     _re_constant(getParam<Real>("re_constant")),
      67          75 :     _well_constant(getParam<Real>("well_constant")),
      68             :     _has_permeability(
      69          74 :         hasMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp") &&
      70          73 :         hasMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")),
      71             :     _has_thermal_conductivity(
      72          66 :         hasMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp") &&
      73          57 :         hasMaterialProperty<std::vector<RealTensorValue>>(
      74             :             "dPorousFlow_thermal_conductivity_qp_dvar")),
      75          25 :     _perm_or_cond(_p_or_t == PorTchoice::pressure
      76          75 :                       ? getMaterialProperty<RealTensorValue>("PorousFlow_permeability_qp")
      77          35 :                       : getMaterialProperty<RealTensorValue>("PorousFlow_thermal_conductivity_qp")),
      78             :     _dperm_or_cond_dvar(
      79          25 :         _p_or_t == PorTchoice::pressure
      80          75 :             ? getMaterialProperty<std::vector<RealTensorValue>>("dPorousFlow_permeability_qp_dvar")
      81             :             : getMaterialProperty<std::vector<RealTensorValue>>(
      82         298 :                   "dPorousFlow_thermal_conductivity_qp_dvar"))
      83             : {
      84          25 :   if (_p_or_t == PorTchoice::pressure && !_has_permeability)
      85             :     mooseError("PorousFlowPeacemanBorehole: You have specified function_of=porepressure, but you "
      86           1 :                "do not have a quadpoint permeability material");
      87          24 :   if (_p_or_t == PorTchoice::temperature && !_has_thermal_conductivity)
      88             :     mooseError("PorousFlowPeacemanBorehole: You have specified function_of=temperature, but you do "
      89           1 :                "not have a quadpoint thermal_conductivity material");
      90             : 
      91             :   // construct the rotation matrix needed to rotate the permeability
      92          23 :   const unsigned int num_pts = _zs.size();
      93          46 :   _rot_matrix.resize(std::max(num_pts - 1, (unsigned)1));
      94          61 :   for (unsigned int i = 0; i + 1 < num_pts; ++i)
      95             :   {
      96         133 :     const RealVectorValue v2(_xs[i + 1] - _xs[i], _ys[i + 1] - _ys[i], _zs[i + 1] - _zs[i]);
      97          19 :     _rot_matrix[i] = RotationMatrix::rotVecToZ(v2);
      98             :   }
      99          23 :   if (num_pts == (unsigned)1)
     100          13 :     _rot_matrix[0] = RotationMatrix::rotVecToZ(_line_direction);
     101          23 : }
     102             : 
     103             : Real
     104      135505 : PorousFlowPeacemanBorehole::wellConstant(const RealTensorValue & perm,
     105             :                                          const RealTensorValue & rot,
     106             :                                          const Real & half_len,
     107             :                                          const Elem * ele,
     108             :                                          const Real & rad) const
     109             : // Peaceman's form for the borehole well constant
     110             : {
     111      135505 :   if (_well_constant > 0)
     112             :     return _well_constant;
     113             : 
     114             :   // rot_perm has its "2" component lying along the half segment.
     115             :   // We want to determine the eigenvectors of rot(0:1, 0:1), since, when
     116             :   // rotated back to the original frame we will determine the element
     117             :   // lengths along these directions
     118      271010 :   const RealTensorValue rot_perm = (rot * perm) * rot.transpose();
     119      135505 :   const Real trace2D = rot_perm(0, 0) + rot_perm(1, 1);
     120      135505 :   const Real det2D = rot_perm(0, 0) * rot_perm(1, 1) - rot_perm(0, 1) * rot_perm(1, 0);
     121      406515 :   const Real sq = std::sqrt(std::max(0.25 * trace2D * trace2D - det2D,
     122      406515 :                                      0.0)); // the std::max accounts for wierdo precision loss
     123      135505 :   const Real eig_val1 = 0.5 * trace2D + sq;
     124      135505 :   const Real eig_val2 = 0.5 * trace2D - sq;
     125             :   RealVectorValue eig_vec1, eig_vec2;
     126      135505 :   if (sq > std::abs(trace2D) * 1E-7) // matrix is not a multiple of the identity (1E-7 accounts for
     127             :                                      // precision in a crude way)
     128             :   {
     129       71344 :     if (rot_perm(1, 0) != 0)
     130             :     {
     131       60368 :       eig_vec1(0) = eig_val1 - rot_perm(1, 1);
     132       60368 :       eig_vec1(1) = rot_perm(1, 0);
     133       60368 :       eig_vec2(0) = eig_val2 - rot_perm(1, 1);
     134       60368 :       eig_vec2(1) = rot_perm(1, 0);
     135             :     }
     136       10976 :     else if (rot_perm(0, 1) != 0)
     137             :     {
     138           0 :       eig_vec1(0) = rot_perm(0, 1);
     139           0 :       eig_vec1(1) = eig_val1 - rot_perm(0, 0);
     140           0 :       eig_vec2(0) = rot_perm(0, 1);
     141           0 :       eig_vec2(1) = eig_val2 - rot_perm(0, 0);
     142             :     }
     143             :     else // off diagonal terms are both zero
     144             :     {
     145       10976 :       eig_vec1(0) = 1.0;
     146       10976 :       eig_vec2(1) = 1.0;
     147             :     }
     148             :   }
     149             :   else // matrix is basically a multiple of the identity
     150             :   {
     151       64161 :     eig_vec1(0) = 1.0;
     152       64161 :     eig_vec2(1) = 1.0;
     153             :   }
     154             : 
     155             :   // finally, rotate these to original frame and normalise
     156      135505 :   eig_vec1 = rot.transpose() * eig_vec1;
     157      135505 :   eig_vec1 /= std::sqrt(eig_vec1 * eig_vec1);
     158      135505 :   eig_vec2 = rot.transpose() * eig_vec2;
     159      135505 :   eig_vec2 /= std::sqrt(eig_vec2 * eig_vec2);
     160             : 
     161             :   // find the "length" of the element in these directions
     162             :   // TODO - maybe better to use variance than max&min
     163             :   Real max1 = eig_vec1 * ele->point(0);
     164             :   Real max2 = eig_vec2 * ele->point(0);
     165             :   Real min1 = max1;
     166             :   Real min2 = max2;
     167             :   Real proj;
     168     2032575 :   for (unsigned int i = 1; i < ele->n_nodes(); i++)
     169             :   {
     170             :     proj = eig_vec1 * ele->point(i);
     171      948535 :     max1 = (max1 < proj) ? proj : max1;
     172      948535 :     min1 = (min1 < proj) ? min1 : proj;
     173             : 
     174             :     proj = eig_vec2 * ele->point(i);
     175      948535 :     max2 = (max2 < proj) ? proj : max2;
     176      948535 :     min2 = (min2 < proj) ? min2 : proj;
     177             :   }
     178      135505 :   const Real ll1 = max1 - min1;
     179      135505 :   const Real ll2 = max2 - min2;
     180             : 
     181             :   Real r0;
     182      135505 :   if (eig_val1 <= 0.0)
     183           0 :     r0 = _re_constant * ll1;
     184      135505 :   else if (eig_val2 <= 0.0)
     185       10976 :     r0 = _re_constant * ll2;
     186             :   else
     187      498116 :     r0 = _re_constant * std::sqrt(std::sqrt(eig_val1 / eig_val2) * std::pow(ll2, 2) +
     188      373587 :                                   std::sqrt(eig_val2 / eig_val1) * std::pow(ll1, 2)) /
     189      124529 :          (std::pow(eig_val1 / eig_val2, 0.25) + std::pow(eig_val2 / eig_val1, 0.25));
     190             : 
     191      135505 :   const Real effective_perm = (det2D >= 0.0 ? std::sqrt(det2D) : 0.0);
     192             : 
     193             :   const Real halfPi = acos(0.0);
     194             : 
     195      135505 :   if (r0 <= rad)
     196             :     mooseError("The effective element size (about 0.2-times-true-ele-size) for an element "
     197             :                "containing a Peaceman-type borehole must be (much) larger than the borehole radius "
     198             :                "for the Peaceman formulation to be correct.  Your element has effective size ",
     199             :                r0,
     200             :                " and the borehole radius is ",
     201             :                rad,
     202           1 :                "\n");
     203             : 
     204      135504 :   return 4 * halfPi * effective_perm * half_len / std::log(r0 / rad);
     205             : }
     206             : 
     207             : Real
     208       49873 : PorousFlowPeacemanBorehole::computeQpBaseOutflow(unsigned current_dirac_ptid) const
     209             : {
     210       99746 :   const Real character = _character.value(_t, _q_point[_qp]);
     211       49873 :   if (character == 0.0)
     212             :     return 0.0;
     213             : 
     214       82946 :   const Real bh_pressure = _p_bot + _unit_weight * (_q_point[_qp] - _bottom_point);
     215       41473 :   const Real pp = ptqp();
     216             : 
     217             :   Real outflow = 0.0; // this is the flow rate from porespace out of the system
     218             : 
     219       41473 :   if (current_dirac_ptid > 0)
     220             :   // contribution from half-segment "behind" this point (must have >1 point for
     221             :   // current_dirac_ptid>0)
     222             :   {
     223       16536 :     if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
     224             :     {
     225             :       // injection, so outflow<0 || production, so outflow>0
     226       33072 :       const Real wc = wellConstant(_perm_or_cond[_qp],
     227             :                                    _rot_matrix[current_dirac_ptid - 1],
     228       16536 :                                    _half_seg_len[current_dirac_ptid - 1],
     229       16536 :                                    _current_elem,
     230       33072 :                                    _rs[current_dirac_ptid]);
     231       16536 :       outflow += wc * (pp - bh_pressure);
     232             :     }
     233             :   }
     234             : 
     235       82946 :   if (current_dirac_ptid + 1 < _zs.size() || _zs.size() == 1)
     236             :   // contribution from half-segment "ahead of" this point, or we only have one point
     237             :   {
     238       33337 :     if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
     239             :     {
     240             :       // injection, so outflow<0 || // production, so outflow>0
     241       66674 :       const Real wc = wellConstant(_perm_or_cond[_qp],
     242             :                                    _rot_matrix[current_dirac_ptid],
     243             :                                    _half_seg_len[current_dirac_ptid],
     244       33337 :                                    _current_elem,
     245       66674 :                                    _rs[current_dirac_ptid]);
     246       33336 :       outflow += wc * (pp - bh_pressure);
     247             :     }
     248             :   }
     249             : 
     250      124416 :   return outflow * _test[_i][_qp] * std::abs(character);
     251             : }
     252             : 
     253             : void
     254       85632 : PorousFlowPeacemanBorehole::computeQpBaseOutflowJacobian(unsigned jvar,
     255             :                                                          unsigned current_dirac_ptid,
     256             :                                                          Real & outflow,
     257             :                                                          Real & outflowp) const
     258             : {
     259       85632 :   outflow = 0.0;
     260       85632 :   outflowp = 0.0;
     261             : 
     262      171264 :   const Real character = _character.value(_t, _q_point[_qp]);
     263       85632 :   if (character == 0.0)
     264             :     return;
     265             : 
     266       77568 :   if (_dictator.notPorousFlowVariable(jvar))
     267             :     return;
     268       77568 :   const unsigned pvar = _dictator.porousFlowVariableNum(jvar);
     269             : 
     270      155136 :   const Real bh_pressure = _p_bot + _unit_weight * (_q_point[_qp] - _bottom_point);
     271       77568 :   const Real pp = ptqp();
     272      155136 :   const Real pp_prime = dptqp(pvar) * _phi[_j][_qp];
     273             : 
     274       77568 :   if (current_dirac_ptid > 0)
     275             :   // contribution from half-segment "behind" this point
     276             :   {
     277       34752 :     if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
     278             :     {
     279             :       // injection, so outflow<0 || // production, so outflow>0
     280       69504 :       const Real wc = wellConstant(_perm_or_cond[_qp],
     281             :                                    _rot_matrix[current_dirac_ptid - 1],
     282       34752 :                                    _half_seg_len[current_dirac_ptid - 1],
     283       34752 :                                    _current_elem,
     284       69504 :                                    _rs[current_dirac_ptid]);
     285       34752 :       outflowp += wc * pp_prime;
     286       34752 :       outflow += wc * (pp - bh_pressure);
     287             :     }
     288             :   }
     289             : 
     290      155136 :   if (current_dirac_ptid < _zs.size() - 1 || _zs.size() == 1)
     291             :   // contribution from half-segment "ahead of" this point
     292             :   {
     293       50880 :     if ((character < 0.0 && pp < bh_pressure) || (character > 0.0 && pp > bh_pressure))
     294             :     {
     295             :       // injection, so outflow<0 || // production, so outflow>0
     296      101760 :       const Real wc = wellConstant(_perm_or_cond[_qp],
     297             :                                    _rot_matrix[current_dirac_ptid],
     298             :                                    _half_seg_len[current_dirac_ptid],
     299       50880 :                                    _current_elem,
     300       50880 :                                    _rs[current_dirac_ptid]);
     301       50880 :       outflowp += wc * pp_prime;
     302       50880 :       outflow += wc * (pp - bh_pressure);
     303             :     }
     304             :   }
     305             : 
     306      232704 :   outflowp *= _test[_i][_qp] * std::abs(character);
     307       77568 :   outflow *= _test[_i][_qp] * std::abs(character);
     308        2499 : }

Generated by: LCOV version 1.11