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 "SamplerBase.h"
11 : #include "IndirectSort.h"
12 : #include "InputParameters.h"
13 : #include "MooseEnum.h"
14 : #include "MooseError.h"
15 : #include "VectorPostprocessor.h"
16 : #include "MooseVariableFieldBase.h"
17 : #include "FEProblemBase.h"
18 : #include "MooseApp.h"
19 : #include "TransientBase.h"
20 :
21 : #include "libmesh/point.h"
22 :
23 : InputParameters
24 31294 : SamplerBase::validParams()
25 : {
26 31294 : InputParameters params = emptyInputParameters();
27 :
28 125176 : params.addRequiredParam<std::string>(
29 : "sort_by",
30 : "What to sort the samples by. Options include 'x', 'y', 'z', 'id', and the name of any of "
31 : "the sampled quantities (which each create a vector of the same name).");
32 :
33 : // The value from this VPP is naturally already on every processor
34 : // TODO: Make this not the case! See #11415
35 31294 : params.set<bool>("_auto_broadcast") = false;
36 :
37 31294 : return params;
38 0 : }
39 :
40 1897 : SamplerBase::SamplerBase(const InputParameters & parameters,
41 : VectorPostprocessor * vpp,
42 1897 : const libMesh::Parallel::Communicator & comm)
43 1897 : : _sampler_params(parameters),
44 1897 : _vpp(vpp),
45 1897 : _comm(comm),
46 1897 : _sampler_transient(dynamic_cast<TransientBase *>(
47 7588 : parameters.getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")
48 1897 : ->getMooseApp()
49 1897 : .getExecutioner())),
50 1897 : _sort_by_index(libMesh::invalid_uint),
51 3794 : _x(vpp->declareVector("x")),
52 3794 : _y(vpp->declareVector("y")),
53 3794 : _z(vpp->declareVector("z")),
54 3794 : _id(vpp->declareVector("id")),
55 5691 : _sort_by(parameters.get<std::string>("sort_by"))
56 : {
57 1897 : if (_sort_by == "x")
58 881 : _sort_by_index = 0;
59 1016 : else if (_sort_by == "y")
60 126 : _sort_by_index = 1;
61 890 : else if (_sort_by == "z")
62 0 : _sort_by_index = 2;
63 890 : else if (_sort_by == "id")
64 877 : _sort_by_index = 3;
65 : // Sort by index for variables determined in setupVariables()
66 1897 : }
67 :
68 : void
69 1879 : SamplerBase::setupVariables(const std::vector<std::string> & variable_names)
70 : {
71 1879 : _variable_names = variable_names;
72 1879 : _values.reserve(variable_names.size());
73 :
74 4509 : for (const auto & variable_name : variable_names)
75 2630 : _values.push_back(&_vpp->declareVector(variable_name));
76 :
77 : // Find the index of the column to sort by for variables
78 1879 : if (_sort_by_index == libMesh::invalid_uint)
79 : {
80 13 : const auto it = std::find(_variable_names.begin(), _variable_names.end(), _sort_by);
81 13 : if (it != _variable_names.end())
82 13 : _sort_by_index = 4 + it - _variable_names.begin();
83 : else
84 0 : mooseError(
85 0 : "The 'sort_by' parameter must be one of x/y/z/id or one of the sampled variable names: " +
86 0 : Moose::stringify(_variable_names));
87 : }
88 1879 : }
89 :
90 : void
91 220695 : SamplerBase::addSample(const Point & p, const Real & id, const std::vector<Real> & values)
92 : {
93 220695 : _x.push_back(p(0));
94 220695 : _y.push_back(p(1));
95 220695 : _z.push_back(p(2));
96 :
97 220695 : _id.push_back(id);
98 :
99 : mooseAssert(values.size() == _variable_names.size(), "Mismatch of variable names to vector size");
100 447907 : for (MooseIndex(values) i = 0; i < values.size(); ++i)
101 227212 : _values[i]->emplace_back(values[i]);
102 :
103 220695 : _curr_num_samples++;
104 220695 : }
105 :
106 : void
107 9430 : SamplerBase::initialize()
108 : {
109 : // Don't reset the vectors if we want to retain history
110 9430 : if (_vpp->containsCompleteHistory() && _comm.rank() == 0)
111 : {
112 : // If we are repeating the timestep due to an aborted solve, we want to throw away the last
113 : // values
114 203 : if (_sampler_transient && !_sampler_transient->lastSolveConverged())
115 : {
116 : // Convenient to allocate a single variable for all the vpp values
117 6 : std::vector<VectorPostprocessorValue *> vec_ptrs = {{&_x, &_y, &_z, &_id}};
118 3 : vec_ptrs.insert(vec_ptrs.end(), _values.begin(), _values.end());
119 : // Erase the elements from the last execution
120 18 : for (auto vec_ptr : vec_ptrs)
121 : {
122 : // Vector may have already been restored, if so, skip the erasure
123 15 : if (_curr_total_samples > vec_ptr->size())
124 : {
125 : mooseAssert(vec_ptr->size() == (_curr_total_samples - _curr_num_samples),
126 : "Number of samples is not what is expected.");
127 1 : continue;
128 : }
129 28 : for (auto ind : _curr_indices)
130 : {
131 : mooseAssert(ind < vec_ptr->size(), "Trying to remove a sample that doesn't exist.");
132 14 : vec_ptr->erase(vec_ptr->begin() + ind);
133 : }
134 : }
135 3 : }
136 203 : _curr_indices.clear();
137 : }
138 : else
139 : {
140 9227 : _x.clear();
141 9227 : _y.clear();
142 9227 : _z.clear();
143 9227 : _id.clear();
144 :
145 9227 : std::for_each(
146 10158 : _values.begin(), _values.end(), [](VectorPostprocessorValue * vec) { vec->clear(); });
147 : }
148 :
149 9430 : _curr_num_samples = 0;
150 9430 : }
151 :
152 : void
153 2224 : SamplerBase::checkForStandardFieldVariableType(const MooseVariableFieldBase * const var_ptr,
154 : const std::string & var_param_name) const
155 : {
156 : // A pointer to a MooseVariableFieldBase should never be SCALAR
157 : mooseAssert(var_ptr->feType().family != SCALAR,
158 : "Scalar variable '" + var_ptr->name() + "' cannot be sampled.");
159 : mooseAssert(dynamic_cast<const MooseObject *>(_vpp), "Should have succeeded");
160 2224 : if (var_ptr->isVector())
161 12 : dynamic_cast<const MooseObject *>(_vpp)->paramError(
162 : var_param_name,
163 : "The variable '",
164 6 : var_ptr->name(),
165 : "' is a vector variable. Sampling those is not currently supported in the "
166 : "framework. It may be supported using a dedicated object in your application. Use "
167 : "'VectorVariableComponentAux' auxkernel to copy those values into a regular field "
168 : "variable");
169 2218 : if (var_ptr->isArray())
170 12 : dynamic_cast<const MooseObject *>(_vpp)->paramError(
171 : var_param_name,
172 : "The variable '",
173 6 : var_ptr->name(),
174 : "' is an array variable. Sampling those is not currently supported in the "
175 : "framework. It may be supported using a dedicated object in your application. Use "
176 : "'ArrayVariableComponent' auxkernel to copy those values into a regular field variable");
177 2212 : }
178 :
179 : void
180 9005 : SamplerBase::finalize()
181 : {
182 : /**
183 : * We have several vectors that all need to be processed in the same way.
184 : * Rather than enumerate them all, let's just create a vector of pointers
185 : * and work on them that way.
186 : */
187 9005 : constexpr auto NUM_ID_VECTORS = 4;
188 :
189 9005 : std::vector<VectorPostprocessorValue *> vec_ptrs;
190 9005 : vec_ptrs.reserve(_values.size() + NUM_ID_VECTORS);
191 : // Initialize the pointer vector with the position and ID vectors
192 18010 : vec_ptrs = {{&_x, &_y, &_z, &_id}};
193 : // Now extend the vector by all the remaining values vector before processing
194 9005 : vec_ptrs.insert(vec_ptrs.end(), _values.begin(), _values.end());
195 :
196 : // Gather up each of the partial vectors
197 54946 : for (auto vec_ptr : vec_ptrs)
198 45941 : _comm.allgather(*vec_ptr, /* identical buffer lengths = */ false);
199 :
200 : // Now create an index vector by using an indirect sort
201 9005 : std::vector<std::size_t> sorted_indices;
202 18010 : Moose::indirectSort(
203 18010 : vec_ptrs[_sort_by_index]->begin(), vec_ptrs[_sort_by_index]->end(), sorted_indices);
204 :
205 : /**
206 : * We now have one sorted vector. The remaining vectors need to be sorted according to that
207 : * vector.
208 : * We'll need a temporary vector to hold values as we map them according to the sorted indices.
209 : * After that, we'll swap the vector contents with the sorted vector to get the values
210 : * back into the original vector.
211 : */
212 : // This vector is used as temp storage to sort each of the remaining vectors according to the
213 : // first
214 9005 : auto vector_length = sorted_indices.size();
215 9005 : VectorPostprocessorValue tmp_vector(vector_length);
216 :
217 : #ifndef NDEBUG
218 : for (const auto vec_ptr : vec_ptrs)
219 : if (vec_ptr->size() != vector_length)
220 : mooseError("Vector length mismatch");
221 : #endif
222 :
223 : // Sort each of the vectors using the same sorted indices
224 54946 : for (auto & vec_ptr : vec_ptrs)
225 : {
226 1605002 : for (MooseIndex(sorted_indices) j = 0; j < sorted_indices.size(); ++j)
227 1559061 : tmp_vector[j] = (*vec_ptr)[sorted_indices[j]];
228 :
229 : // Swap vector storage with sorted vector
230 45941 : vec_ptr->swap(tmp_vector);
231 : }
232 :
233 : // Gather the indices of samples from the last execution
234 : // Used to determine which parts of the vector need to be erased if a solve fails
235 9005 : if (_vpp->containsCompleteHistory())
236 : {
237 278 : _comm.sum(_curr_num_samples);
238 278 : if (_comm.rank() == 0)
239 : {
240 203 : _curr_indices.insert(sorted_indices.end() - _curr_num_samples, sorted_indices.end());
241 203 : _curr_total_samples = vec_ptrs[0]->size();
242 : }
243 : }
244 9005 : }
245 :
246 : void
247 422 : SamplerBase::threadJoin(const SamplerBase & y)
248 : {
249 422 : _x.insert(_x.end(), y._x.begin(), y._x.end());
250 422 : _y.insert(_y.end(), y._y.begin(), y._y.end());
251 422 : _z.insert(_z.end(), y._z.begin(), y._z.end());
252 :
253 422 : _id.insert(_id.end(), y._id.begin(), y._id.end());
254 :
255 856 : for (MooseIndex(_variable_names) i = 0; i < _variable_names.size(); i++)
256 434 : _values[i]->insert(_values[i]->end(), y._values[i]->begin(), y._values[i]->end());
257 :
258 422 : _curr_num_samples += y._curr_num_samples;
259 422 : }
|
