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