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 "PorousFlowLineGeometry.h"
11 : #include "RayTracing.h"
12 : #include "LineSegment.h"
13 : #include "libmesh/utility.h"
14 :
15 : #include <fstream>
16 :
17 : InputParameters
18 2014 : PorousFlowLineGeometry::validParams()
19 : {
20 2014 : InputParameters params = DiracKernel::validParams();
21 4028 : 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 4028 : 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 4028 : 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 4028 : 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 4028 : 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 6042 : params.addRangeCheckedParam<Real>(
52 : "line_length",
53 4028 : 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 2014 : params.addParam<RealVectorValue>(
57 : "line_direction",
58 2014 : 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 4028 : 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 2014 : 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 2014 : return params;
67 0 : }
68 :
69 1131 : PorousFlowLineGeometry::PorousFlowLineGeometry(const InputParameters & parameters)
70 : : DiracKernel(parameters),
71 : ReporterInterface(this),
72 1131 : _line_length(getParam<Real>("line_length")),
73 2262 : _line_direction(getParam<RealVectorValue>("line_direction")),
74 3393 : _point_file(getParam<std::string>("point_file")),
75 2412 : _x_coord(isParamValid("x_coord_reporter") ? &getReporterValue<std::vector<Real>>(
76 : "x_coord_reporter", REPORTER_MODE_REPLICATED)
77 : : &_xs),
78 2412 : _y_coord(isParamValid("y_coord_reporter") ? &getReporterValue<std::vector<Real>>(
79 : "y_coord_reporter", REPORTER_MODE_REPLICATED)
80 : : &_ys),
81 2412 : _z_coord(isParamValid("z_coord_reporter") ? &getReporterValue<std::vector<Real>>(
82 : "z_coord_reporter", REPORTER_MODE_REPLICATED)
83 : : &_zs),
84 1131 : _weight(isParamValid("weight_reporter")
85 1281 : ? &getReporterValue<std::vector<Real>>("weight_reporter", REPORTER_MODE_REPLICATED)
86 : : &_rs),
87 3393 : _usingReporter(isParamValid("x_coord_reporter"))
88 : {
89 1131 : statefulPropertiesAllowed(true);
90 :
91 : const int checkInputFormat =
92 2262 : int(isParamValid("line_base")) + int(!_point_file.empty()) + int(_usingReporter);
93 :
94 1131 : if (checkInputFormat > 1)
95 0 : paramError("point_file",
96 : "PorousFlowLineGeometry: must specify only one of 'point_file' or 'line_base' or "
97 : "reporter based input");
98 1131 : else if (checkInputFormat == 0)
99 0 : paramError(
100 : "point_file",
101 : "PorousFlowLineGeometry: must specify at least one of 'point_file' or 'line_base' or "
102 : "reporter based input");
103 1131 : }
104 :
105 : void
106 1099 : PorousFlowLineGeometry::initialSetup()
107 : {
108 1099 : if (!_point_file.empty())
109 : {
110 : // open file
111 1025 : std::ifstream file(_point_file.c_str());
112 1025 : if (!file.good())
113 2 : 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 2993 : while (parseNextLineReals(file, scratch))
118 : {
119 1970 : if (scratch.size() >= 2)
120 : {
121 1767 : _rs.push_back(scratch[0]);
122 1767 : _xs.push_back(scratch[1]);
123 1767 : if (scratch.size() >= 3)
124 1767 : _ys.push_back(scratch[2]);
125 : else
126 0 : _ys.push_back(0.0);
127 1767 : if (scratch.size() >= 4)
128 1213 : _zs.push_back(scratch[3]);
129 : else
130 554 : _zs.push_back(0.0);
131 : }
132 : }
133 1023 : file.close();
134 1023 : calcLineLengths();
135 1021 : }
136 74 : else if (_usingReporter)
137 : {
138 48 : if (_weight->size() != _x_coord->size() || _weight->size() != _y_coord->size() ||
139 46 : _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 2 : "weight size = " +
148 4 : std::to_string(_weight->size()) +
149 6 : "; x_coord size = " + std::to_string(_x_coord->size()) +
150 6 : "; y_coord size = " + std::to_string(_y_coord->size()) +
151 4 : "; z_coord size = " + std::to_string(_z_coord->size());
152 :
153 2 : mooseError(errMsg);
154 : }
155 :
156 134 : for (std::size_t i = 0; i < _x_coord->size(); ++i)
157 : {
158 88 : _rs.push_back(_weight->at(i));
159 88 : _xs.push_back(_x_coord->at(i));
160 88 : _ys.push_back(_y_coord->at(i));
161 88 : _zs.push_back(_z_coord->at(i));
162 : }
163 46 : calcLineLengths();
164 : }
165 : else
166 : {
167 78 : _line_base = getParam<std::vector<Real>>("line_base");
168 26 : if (_line_base.size() != _mesh.dimension() + 1)
169 0 : paramError("line_base",
170 : "PorousFlowLineGeometry: wrong number of arguments - got ",
171 : _line_base.size(),
172 : ", expected ",
173 0 : _mesh.dimension() + 1,
174 : " '<weight> <x> [<y> [z]]'");
175 :
176 40 : for (size_t i = _line_base.size(); i < 4; i++)
177 14 : _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 26 : Point start(_line_base[1], _line_base[2], _line_base[3]);
181 26 : Point end = start + _line_length * _line_direction / _line_direction.norm();
182 26 : auto pl = _subproblem.mesh().getPointLocator();
183 26 : pl->enable_out_of_mesh_mode();
184 :
185 26 : auto * elem = (*pl)(start);
186 26 : auto elem_id = elem ? elem->id() : libMesh::DofObject::invalid_id;
187 26 : _communicator.min(elem_id);
188 26 : if (elem_id == libMesh::DofObject::invalid_id)
189 0 : paramError("line_base", "base point ", start, " lies outside the mesh");
190 :
191 26 : elem = (*pl)(end);
192 26 : elem_id = elem ? elem->id() : libMesh::DofObject::invalid_id;
193 26 : _communicator.min(elem_id);
194 26 : if (elem_id == libMesh::DofObject::invalid_id)
195 2 : paramError("line_length", "length causes end point ", end, " to lie outside the mesh");
196 :
197 24 : regenPoints();
198 24 : }
199 1089 : }
200 :
201 : void
202 1120 : PorousFlowLineGeometry::calcLineLengths()
203 : {
204 1120 : const int num_pts = _zs.size();
205 1120 : if (num_pts == 0)
206 2 : mooseError("PorousFlowLineGeometry: No points found in input.\nIf using reporters, make sure "
207 : "they have data.");
208 1118 : _bottom_point(0) = _xs[num_pts - 1];
209 1118 : _bottom_point(1) = _ys[num_pts - 1];
210 1118 : _bottom_point(2) = _zs[num_pts - 1];
211 :
212 : // construct the line-segment lengths between each point
213 1118 : _half_seg_len.clear();
214 1781 : _half_seg_len.resize(std::max(num_pts - 1, 1));
215 1943 : for (unsigned int i = 0; i + 1 < _xs.size(); ++i)
216 : {
217 827 : _half_seg_len[i] = 0.5 * std::sqrt(Utility::pow<2>(_xs[i + 1] - _xs[i]) +
218 827 : Utility::pow<2>(_ys[i + 1] - _ys[i]) +
219 827 : Utility::pow<2>(_zs[i + 1] - _zs[i]));
220 827 : if (_half_seg_len[i] == 0)
221 2 : 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 1116 : if (num_pts == 1)
230 663 : _half_seg_len[0] = _line_length;
231 1116 : }
232 :
233 : void
234 51 : PorousFlowLineGeometry::regenPoints()
235 : {
236 51 : if (!_point_file.empty() || _usingReporter)
237 0 : return;
238 :
239 : // recalculate the auto-generated points:
240 51 : _rs.clear();
241 51 : _xs.clear();
242 51 : _ys.clear();
243 51 : _zs.clear();
244 :
245 51 : Point p0(_line_base[1], _line_base[2], _line_base[3]);
246 51 : Point p1 = p0 + _line_length * _line_direction / _line_direction.norm();
247 :
248 : // add point for each cell the line passes through
249 51 : auto ploc = _mesh.getPointLocator();
250 : std::vector<Elem *> elems;
251 : std::vector<LineSegment> segs;
252 51 : Moose::elementsIntersectedByLine(p0, p1, _mesh, *ploc, elems, segs);
253 141 : 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 0 : if (seg.start() == seg.end())
258 0 : continue;
259 :
260 : auto middle = (seg.start() + seg.end()) * 0.5;
261 90 : _rs.push_back(_line_base[0]);
262 90 : _xs.push_back(middle(0));
263 90 : _ys.push_back(middle(1));
264 90 : _zs.push_back(middle(2));
265 : }
266 :
267 : // make the start point be the line base point
268 51 : _rs.front() = _line_base[0];
269 51 : _xs.front() = p0(0);
270 51 : _ys.front() = p0(1);
271 51 : _zs.front() = p0(2);
272 :
273 : // force the end point only if our line traverses more than one element
274 51 : if (segs.size() > 1)
275 : {
276 30 : _rs.back() = _line_base[0];
277 30 : _xs.back() = p1(0);
278 30 : _ys.back() = p1(1);
279 30 : _zs.back() = p1(2);
280 : }
281 51 : calcLineLengths();
282 51 : }
283 :
284 : void
285 27 : PorousFlowLineGeometry::meshChanged()
286 : {
287 27 : DiracKernel::meshChanged();
288 27 : regenPoints();
289 27 : }
290 :
291 : bool
292 2993 : 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 2993 : myvec.clear();
297 : bool gotline(false);
298 2993 : if (getline(ifs, line))
299 : {
300 : gotline = true;
301 :
302 : // Harvest floats separated by whitespace
303 1970 : std::istringstream iss(line);
304 : Real f;
305 8484 : while (iss >> f)
306 : {
307 6514 : myvec.push_back(f);
308 : }
309 1970 : }
310 2993 : return gotline;
311 : }
312 :
313 : void
314 62179 : PorousFlowLineGeometry::addPoints()
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 182532 : for (unsigned int i = 0; i < _x_coord->size(); i++)
320 120353 : addPoint(Point(_x_coord->at(i), _y_coord->at(i), _z_coord->at(i)), i);
321 62179 : }
|
