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 "RichardsPolyLineSink.h"
      11             : 
      12             : #include <fstream>
      13             : 
      14             : registerMooseObject("RichardsApp", RichardsPolyLineSink);
      15             : 
      16             : InputParameters
      17          27 : RichardsPolyLineSink::validParams()
      18             : {
      19          27 :   InputParameters params = DiracKernel::validParams();
      20          54 :   params.addRequiredParam<std::vector<Real>>(
      21             :       "pressures", "Tuple of pressure values.  Must be monotonically increasing.");
      22          54 :   params.addRequiredParam<std::vector<Real>>(
      23             :       "fluxes",
      24             :       "Tuple of flux values (measured in kg.m^-3.s^-1).  A piecewise-linear fit is "
      25             :       "performed to the (pressure,flux) pairs to obtain the flux at any arbitrary "
      26             :       "pressure.  If a quad-point pressure is less than the first pressure value, the "
      27             :       "first flux value is used.  If quad-point pressure exceeds the final pressure "
      28             :       "value, the final flux value is used.  This flux is OUT of the medium: hence "
      29             :       "positive values of flux means this will be a SINK, while negative values indicate "
      30             :       "this flux will be a SOURCE.");
      31          54 :   params.addRequiredParam<FileName>(
      32             :       "point_file",
      33             :       "The file containing the coordinates of the point sinks that will approximate "
      34             :       "the polyline.  Each line in the file must contain a space-separated "
      35             :       "coordinate.  Note that you will get segementation faults if your points do "
      36             :       "not lie within your mesh!");
      37          54 :   params.addRequiredParam<UserObjectName>(
      38             :       "SumQuantityUO",
      39             :       "User Object of type=RichardsSumQuantity in which to place the total "
      40             :       "outflow from the polylinesink for each time step.");
      41          54 :   params.addRequiredParam<UserObjectName>(
      42             :       "richardsVarNames_UO", "The UserObject that holds the list of Richards variable names.");
      43          27 :   params.addClassDescription("Approximates a polyline sink in the mesh by using a number of point "
      44             :                              "sinks whose positions are read from a file");
      45          27 :   return params;
      46           0 : }
      47             : 
      48          14 : RichardsPolyLineSink::RichardsPolyLineSink(const InputParameters & parameters)
      49             :   : DiracKernel(parameters),
      50          14 :     _total_outflow_mass(
      51          14 :         const_cast<RichardsSumQuantity &>(getUserObject<RichardsSumQuantity>("SumQuantityUO"))),
      52          56 :     _sink_func(getParam<std::vector<Real>>("pressures"), getParam<std::vector<Real>>("fluxes")),
      53          28 :     _point_file(getParam<FileName>("point_file")),
      54          14 :     _richards_name_UO(getUserObject<RichardsVarNames>("richardsVarNames_UO")),
      55          14 :     _pvar(_richards_name_UO.richards_var_num(_var.number())),
      56          28 :     _pp(getMaterialProperty<std::vector<Real>>("porepressure")),
      57          42 :     _dpp_dv(getMaterialProperty<std::vector<std::vector<Real>>>("dporepressure_dv"))
      58             : {
      59             :   // open file
      60          14 :   std::ifstream file(_point_file.c_str());
      61          14 :   if (!file.good())
      62           2 :     mooseError("Error opening file '" + _point_file + "' from RichardsPolyLineSink.");
      63             : 
      64             :   std::vector<Real> scratch;
      65          27 :   while (parseNextLineReals(file, scratch))
      66             :   {
      67          14 :     if (scratch.size() >= 1)
      68             :     {
      69          14 :       _xs.push_back(scratch[0]);
      70          14 :       if (scratch.size() >= 2)
      71          12 :         _ys.push_back(scratch[1]);
      72             :       else
      73           2 :         _ys.push_back(0.0);
      74             : 
      75          14 :       if (scratch.size() >= 3)
      76          10 :         _zs.push_back(scratch[2]);
      77             :       else
      78           4 :         _zs.push_back(0.0);
      79             :     }
      80             :   }
      81             : 
      82          13 :   file.close();
      83             : 
      84             :   // To correctly compute the Jacobian terms,
      85             :   // tell MOOSE that this DiracKernel depends on all the Richards Vars
      86          13 :   const std::vector<MooseVariableFEBase *> & coupled_vars = _richards_name_UO.getCoupledMooseVars();
      87          28 :   for (unsigned int i = 0; i < coupled_vars.size(); i++)
      88          30 :     addMooseVariableDependency(coupled_vars[i]);
      89          13 : }
      90             : 
      91             : bool
      92          27 : RichardsPolyLineSink::parseNextLineReals(std::ifstream & ifs, std::vector<Real> & myvec)
      93             : // reads a space-separated line of floats from ifs and puts in myvec
      94             : {
      95             :   std::string line;
      96          27 :   myvec.clear();
      97             :   bool gotline(false);
      98          27 :   if (getline(ifs, line))
      99             :   {
     100             :     gotline = true;
     101             : 
     102             :     // Harvest floats separated by whitespace
     103          14 :     std::istringstream iss(line);
     104             :     Real f;
     105          50 :     while (iss >> f)
     106             :     {
     107          36 :       myvec.push_back(f);
     108             :     }
     109          14 :   }
     110          27 :   return gotline;
     111             : }
     112             : 
     113             : void
     114        2162 : RichardsPolyLineSink::addPoints()
     115             : {
     116        2162 :   _total_outflow_mass.zero();
     117             : 
     118             :   // Add point using the unique ID "i", let the DiracKernel take
     119             :   // care of the caching.  This should be fast after the first call,
     120             :   // as long as the points don't move around.
     121        4337 :   for (unsigned int i = 0; i < _zs.size(); i++)
     122        2175 :     addPoint(Point(_xs[i], _ys[i], _zs[i]), i);
     123        2162 : }
     124             : 
     125             : Real
     126       13448 : RichardsPolyLineSink::computeQpResidual()
     127             : {
     128       13448 :   Real test_fcn = _test[_i][_qp];
     129       13448 :   Real flow = test_fcn * _sink_func.sample(_pp[_qp][_pvar]);
     130       13448 :   _total_outflow_mass.add(flow * _dt);
     131       13448 :   return flow;
     132             : }
     133             : 
     134             : Real
     135       31616 : RichardsPolyLineSink::computeQpJacobian()
     136             : {
     137       31616 :   Real test_fcn = _test[_i][_qp];
     138       31616 :   return test_fcn * _sink_func.sampleDerivative(_pp[_qp][_pvar]) * _dpp_dv[_qp][_pvar][_pvar] *
     139       31616 :          _phi[_j][_qp];
     140             : }
     141             : 
     142             : Real
     143        3200 : RichardsPolyLineSink::computeQpOffDiagJacobian(unsigned int jvar)
     144             : {
     145        3200 :   if (_richards_name_UO.not_richards_var(jvar))
     146             :     return 0.0;
     147        3200 :   unsigned int dvar = _richards_name_UO.richards_var_num(jvar);
     148        3200 :   Real test_fcn = _test[_i][_qp];
     149        3200 :   return test_fcn * _sink_func.sampleDerivative(_pp[_qp][_pvar]) * _dpp_dv[_qp][_pvar][dvar] *
     150        3200 :          _phi[_j][_qp];
     151             : }