www.mooseframework.org
PorousFlowLineGeometry.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "PorousFlowLineGeometry.h"
11 #include "RayTracing.h"
12 #include "LineSegment.h"
13 #include "libmesh/utility.h"
14 
15 #include <fstream>
16 
19 {
21  params.addParam<std::string>(
22  "point_file",
23  "",
24  "The file containing the coordinates of the points and their weightings that approximate the "
25  "line sink. The physical meaning of the weightings depend on the scenario, eg, they may be "
26  "borehole radii. Each line in the file must contain a space-separated weight and "
27  "coordinate, viz r x y z. For boreholes, the last point in the file is defined as the "
28  "borehole bottom, where the borehole pressure is bottom_pressure. If your file contains "
29  "just one point, you must also specify the line_length and line_direction parameters. Note "
30  "that you will get segementation faults if your points do not lie within your mesh!");
31  params.addParam<ReporterName>(
32  "x_coord_reporter",
33  "reporter x-coordinate name of line sink. This uses the reporter syntax <reporter>/<name>. "
34  "Each point must adhere to the same requirements as those that would be given if using "
35  "point_file ");
36  params.addParam<ReporterName>(
37  "y_coord_reporter",
38  "reporter y-coordinate name of line sink. This uses the reporter syntax <reporter>/<name>. "
39  "Each point must adhere to the same requirements as those that would be given if using "
40  "point_file ");
41  params.addParam<ReporterName>(
42  "z_coord_reporter",
43  "reporter z-coordinate name of line sink. This uses the reporter syntax <reporter>/<name>. "
44  "Each point must adhere to the same requirements as those that would be given if using "
45  "point_file ");
46  params.addParam<ReporterName>(
47  "weight_reporter",
48  "reporter weight name of line sink. This uses the reporter syntax <reporter>/<name>. "
49  "Each point must adhere to the same requirements as those that would be given if using "
50  "point_file ");
51  params.addRangeCheckedParam<Real>(
52  "line_length",
53  0.0,
54  "line_length>=0",
55  "Line length. Note this is only used if there is only one point in the point_file.");
56  params.addParam<RealVectorValue>(
57  "line_direction",
58  RealVectorValue(0.0, 0.0, 1.0),
59  "Line direction. Note this is only used if there is only one point in the point_file.");
60  params.addParam<std::vector<Real>>(
61  "line_base",
62  "Line base point x,y,z coordinates. This is the same format as a single-line point_file. "
63  "Note this is only used if there is no point file specified.");
64  params.addClassDescription("Approximates a polyline sink in the mesh using a number of Dirac "
65  "point sinks with given weightings that are read from a file");
66  return params;
67 }
68 
70  : DiracKernel(parameters),
71  ReporterInterface(this),
72  _line_length(getParam<Real>("line_length")),
73  _line_direction(getParam<RealVectorValue>("line_direction")),
74  _point_file(getParam<std::string>("point_file")),
75  _x_coord(isParamValid("x_coord_reporter") ? &getReporterValue<std::vector<Real>>(
76  "x_coord_reporter", REPORTER_MODE_REPLICATED)
77  : &_xs),
78  _y_coord(isParamValid("y_coord_reporter") ? &getReporterValue<std::vector<Real>>(
79  "y_coord_reporter", REPORTER_MODE_REPLICATED)
80  : &_ys),
81  _z_coord(isParamValid("z_coord_reporter") ? &getReporterValue<std::vector<Real>>(
82  "z_coord_reporter", REPORTER_MODE_REPLICATED)
83  : &_zs),
84  _weight(isParamValid("weight_reporter")
85  ? &getReporterValue<std::vector<Real>>("weight_reporter", REPORTER_MODE_REPLICATED)
86  : &_rs),
87  _usingReporter(isParamValid("x_coord_reporter"))
88 {
90 
91  const int checkInputFormat =
92  int(isParamValid("line_base")) + int(!_point_file.empty()) + int(_usingReporter);
93 
94  if (checkInputFormat > 1)
95  paramError("point_file",
96  "PorousFlowLineGeometry: must specify only one of 'point_file' or 'line_base' or "
97  "reporter based input");
98  else if (checkInputFormat == 0)
99  paramError(
100  "point_file",
101  "PorousFlowLineGeometry: must specify at least one of 'point_file' or 'line_base' or "
102  "reporter based input");
103 }
104 
105 void
107 {
108  if (!_point_file.empty())
109  {
110  // open file
111  std::ifstream file(_point_file.c_str());
112  if (!file.good())
113  paramError("point_file", "PorousFlowLineGeometry: Error opening file " + _point_file);
114 
115  // construct the arrays of weight, x, y and z
116  std::vector<Real> scratch;
117  while (parseNextLineReals(file, scratch))
118  {
119  if (scratch.size() >= 2)
120  {
121  _rs.push_back(scratch[0]);
122  _xs.push_back(scratch[1]);
123  if (scratch.size() >= 3)
124  _ys.push_back(scratch[2]);
125  else
126  _ys.push_back(0.0);
127  if (scratch.size() >= 4)
128  _zs.push_back(scratch[3]);
129  else
130  _zs.push_back(0.0);
131  }
132  }
133  file.close();
134  calcLineLengths();
135  }
136  else if (_usingReporter)
137  {
138  if (_weight->size() != _x_coord->size() || _weight->size() != _y_coord->size() ||
139  _weight->size() != _z_coord->size())
140  {
141  std::string errMsg =
142  "The value and coordinate vectors are a different size. \n"
143  "There must be one value per coordinate. If the sizes are \n"
144  "zero, the reporter or reporter may not have been initialized with data \n"
145  "before the Dirac Kernel is called. \n"
146  "Try setting \"execute_on = timestep_begin\" in the reporter being read. \n"
147  "weight size = " +
148  std::to_string(_weight->size()) +
149  "; x_coord size = " + std::to_string(_x_coord->size()) +
150  "; y_coord size = " + std::to_string(_y_coord->size()) +
151  "; z_coord size = " + std::to_string(_z_coord->size());
152 
153  mooseError(errMsg);
154  }
155 
156  for (std::size_t i = 0; i < _x_coord->size(); ++i)
157  {
158  _rs.push_back(_weight->at(i));
159  _xs.push_back(_x_coord->at(i));
160  _ys.push_back(_y_coord->at(i));
161  _zs.push_back(_z_coord->at(i));
162  }
163  calcLineLengths();
164  }
165  else
166  {
167  _line_base = getParam<std::vector<Real>>("line_base");
168  if (_line_base.size() != _mesh.dimension() + 1)
169  paramError("line_base",
170  "PorousFlowLineGeometry: wrong number of arguments - got ",
171  _line_base.size(),
172  ", expected ",
173  _mesh.dimension() + 1,
174  " '<weight> <x> [<y> [z]]'");
175 
176  for (size_t i = _line_base.size(); i < 4; i++)
177  _line_base.push_back(0); // fill out zeros up to weight+3 dimensions
178 
179  // make sure line base point is inside the mesh
180  Point start(_line_base[1], _line_base[2], _line_base[3]);
181  Point end = start + _line_length * _line_direction / _line_direction.norm();
182  auto pl = _subproblem.mesh().getPointLocator();
183  pl->enable_out_of_mesh_mode();
184 
185  auto * elem = (*pl)(start);
186  auto elem_id = elem ? elem->id() : DofObject::invalid_id;
187  _communicator.min(elem_id);
188  if (elem_id == DofObject::invalid_id)
189  paramError("line_base", "base point ", start, " lies outside the mesh");
190 
191  elem = (*pl)(end);
192  elem_id = elem ? elem->id() : DofObject::invalid_id;
193  _communicator.min(elem_id);
194  if (elem_id == DofObject::invalid_id)
195  paramError("line_length", "length causes end point ", end, " to lie outside the mesh");
196 
197  regenPoints();
198  }
199 }
200 
201 void
203 {
204  const int num_pts = _zs.size();
205  if (num_pts == 0)
206  mooseError("PorousFlowLineGeometry: No points found in input.\nIf using reporters, make sure "
207  "they have data.");
208  _bottom_point(0) = _xs[num_pts - 1];
209  _bottom_point(1) = _ys[num_pts - 1];
210  _bottom_point(2) = _zs[num_pts - 1];
211 
212  // construct the line-segment lengths between each point
213  _half_seg_len.clear();
214  _half_seg_len.resize(std::max(num_pts - 1, 1));
215  for (unsigned int i = 0; i + 1 < _xs.size(); ++i)
216  {
217  _half_seg_len[i] = 0.5 * std::sqrt(Utility::pow<2>(_xs[i + 1] - _xs[i]) +
218  Utility::pow<2>(_ys[i + 1] - _ys[i]) +
219  Utility::pow<2>(_zs[i + 1] - _zs[i]));
220  if (_half_seg_len[i] == 0)
221  mooseError("PorousFlowLineGeometry: zero-segment length detected at (x,y,z) = ",
222  _xs[i],
223  " ",
224  _ys[i],
225  " ",
226  _zs[i],
227  "\n");
228  }
229  if (num_pts == 1)
231 }
232 
233 void
235 {
236  if (!_point_file.empty() || _usingReporter)
237  return;
238 
239  // recalculate the auto-generated points:
240  _rs.clear();
241  _xs.clear();
242  _ys.clear();
243  _zs.clear();
244 
245  Point p0(_line_base[1], _line_base[2], _line_base[3]);
246  Point p1 = p0 + _line_length * _line_direction / _line_direction.norm();
247 
248  // add point for each cell the line passes through
249  auto ploc = _mesh.getPointLocator();
250  std::vector<Elem *> elems;
251  std::vector<LineSegment> segs;
252  Moose::elementsIntersectedByLine(p0, p1, _mesh, *ploc, elems, segs);
253  for (size_t i = 0; i < segs.size(); i++)
254  {
255  // elementsIntersectedByLine sometimes returns segments with coincident points - check for this:
256  auto & seg = segs[i];
257  if (seg.start() == seg.end())
258  continue;
259 
260  auto middle = (seg.start() + seg.end()) * 0.5;
261  _rs.push_back(_line_base[0]);
262  _xs.push_back(middle(0));
263  _ys.push_back(middle(1));
264  _zs.push_back(middle(2));
265  }
266 
267  // make the start point be the line base point
268  _rs.front() = _line_base[0];
269  _xs.front() = p0(0);
270  _ys.front() = p0(1);
271  _zs.front() = p0(2);
272 
273  // force the end point only if our line traverses more than one element
274  if (segs.size() > 1)
275  {
276  _rs.back() = _line_base[0];
277  _xs.back() = p1(0);
278  _ys.back() = p1(1);
279  _zs.back() = p1(2);
280  }
281  calcLineLengths();
282 }
283 
284 void
286 {
288  regenPoints();
289 }
290 
291 bool
292 PorousFlowLineGeometry::parseNextLineReals(std::ifstream & ifs, std::vector<Real> & myvec)
293 // reads a space-separated line of floats from ifs and puts in myvec
294 {
295  std::string line;
296  myvec.clear();
297  bool gotline(false);
298  if (getline(ifs, line))
299  {
300  gotline = true;
301 
302  // Harvest floats separated by whitespace
303  std::istringstream iss(line);
304  Real f;
305  while (iss >> f)
306  {
307  myvec.push_back(f);
308  }
309  }
310  return gotline;
311 }
312 
313 void
315 {
316  // Add point using the unique ID "i", let the DiracKernel take
317  // care of the caching. This should be fast after the first call,
318  // as long as the points don't move around.
319  for (unsigned int i = 0; i < _x_coord->size(); i++)
320  addPoint(Point(_x_coord->at(i), _y_coord->at(i), _z_coord->at(i)), i);
321 }
MooseMesh & _mesh
virtual MooseMesh & mesh()=0
virtual void addPoints() override
Add Dirac Points to the line sink.
const std::vector< Real > *const _weight
auto norm() const -> decltype(std::norm(Real()))
const std::vector< Real > *const _y_coord
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
virtual void meshChanged() override
regenerate points in each cell if using line_base
bool parseNextLineReals(std::ifstream &ifs, std::vector< Real > &myvec)
Reads a space-separated line of floats from ifs and puts in myvec.
virtual void meshChanged() override
const Real _line_length
Line length. This is only used if there is only one borehole point.
void addPoint(const Elem *elem, Point p, unsigned id=libMesh::invalid_uint)
std::vector< Real > _ys
y points of the borehole
const Parallel::Communicator & _communicator
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
bool isParamValid(const std::string &name) const
std::vector< Real > _half_seg_len
0.5*(length of polyline segments between points)
const RealVectorValue _line_direction
Line direction. This is only used if there is only one borehole point.
void min(const T &r, T &o, Request &req) const
PorousFlowLineGeometry(const InputParameters &parameters)
virtual unsigned int dimension() const
Real f(Real x)
Test function for Brents method.
virtual void initialSetup() override
std::vector< Real > _zs
z points of borehole
SubProblem & _subproblem
void paramError(const std::string &param, Args... args) const
std::vector< Real > _line_base
alternative (to the point file data) line weight and start point.
const std::vector< Real > *const _x_coord
static InputParameters validParams()
Creates a new PorousFlowLineGeometry This reads the file containing the lines of the form weight x y ...
const std::vector< Real > *const _z_coord
const std::string _point_file
File defining the geometry of the borehole.
void statefulPropertiesAllowed(bool)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void mooseError(Args &&... args) const
void addClassDescription(const std::string &doc_string)
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
const ReporterMode REPORTER_MODE_REPLICATED
std::vector< Real > _xs
x points of the borehole
std::vector< Real > _rs
Radii of the borehole.
Point _bottom_point
The bottom point of the borehole (where bottom_pressure is defined)
virtual std::unique_ptr< PointLocatorBase > getPointLocator() const
void ErrorVector unsigned int
void elementsIntersectedByLine(const Point &p0, const Point &p1, const MeshBase &mesh, const PointLocatorBase &point_locator, std::vector< Elem * > &intersected_elems, std::vector< LineSegment > &segments)
static InputParameters validParams()