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 "NearestReporterCoordinatesFunction.h"
11 :
12 : registerMooseObject("OptimizationApp", NearestReporterCoordinatesFunction);
13 :
14 : InputParameters
15 456 : NearestReporterCoordinatesFunction::validParams()
16 : {
17 456 : InputParameters params = OptimizationFunction::validParams();
18 456 : params.addClassDescription(
19 : "This Function finds the nearest point in the specified vectors of coordinates and returns "
20 : "the values specified in the vector of values at the index of the nearest point. All the "
21 : "vectors must be specified using either vector postprocessors or reporter vectors. This "
22 : "function interpolates linearly in time with transient data.");
23 912 : params.addParam<ReporterName>("x_coord_name",
24 : "Name of vector-postprocessor or reporter vector containing "
25 : "x-coordinate of points, default is assumed to be all 0s.");
26 912 : params.addParam<ReporterName>("y_coord_name",
27 : "Name of vector-postprocessor or reporter vector containing "
28 : "y-coordinate of points, default is assumed to be all 0s.");
29 912 : params.addParam<ReporterName>("z_coord_name",
30 : "Name of vector-postprocessor or reporter vector containing "
31 : "z-coordinate of points, default is assumed to be all 0s.");
32 912 : params.addParam<ReporterName>("time_name",
33 : "Name of vector-postprocessor or reporter vector containing time, "
34 : "default is assumed to be all 0s.");
35 912 : params.addRequiredParam<ReporterName>(
36 : "value_name", "Name of vector-postprocessor or reporter vector containing value data.");
37 456 : return params;
38 0 : }
39 :
40 238 : NearestReporterCoordinatesFunction::NearestReporterCoordinatesFunction(
41 238 : const InputParameters & parameters)
42 : : OptimizationFunction(parameters),
43 : ReporterInterface(this),
44 476 : _coordx(isParamValid("x_coord_name") ? getReporterValue<std::vector<Real>>("x_coord_name")
45 : : _empty_vec),
46 659 : _coordy(isParamValid("y_coord_name") ? getReporterValue<std::vector<Real>>("y_coord_name")
47 : : _empty_vec),
48 535 : _coordz(isParamValid("z_coord_name") ? getReporterValue<std::vector<Real>>("z_coord_name")
49 : : _empty_vec),
50 513 : _coordt(isParamValid("time_name") ? getReporterValue<std::vector<Real>>("time_name")
51 : : _empty_vec),
52 476 : _values(getReporterValue<std::vector<Real>>("value_name"))
53 : {
54 238 : }
55 :
56 : Real
57 46112079 : NearestReporterCoordinatesFunction::value(Real t, const Point & p) const
58 : {
59 46112079 : const std::array<std::pair<Real, std::size_t>, 2> tv = findNearestPoint(t, p);
60 :
61 : // If the indices are equal then t is either less than the minimum time in the data
62 : // or greater than the maximum time. In which case will extrapolate using a constant
63 : // value.
64 46112064 : if (tv[0].second == tv[1].second)
65 46106880 : return _values[tv[0].second];
66 :
67 5184 : const Real told = tv[0].first;
68 5184 : const Real tnew = tv[1].first;
69 5184 : const Real vold = _values[tv[0].second];
70 5184 : const Real vnew = _values[tv[1].second];
71 5184 : return vold + (vnew - vold) * (t - told) / (tnew - told);
72 : }
73 :
74 : RealGradient
75 0 : NearestReporterCoordinatesFunction::gradient(Real /*t*/, const Point & /*p*/) const
76 : {
77 0 : return RealGradient(0, 0, 0);
78 : }
79 :
80 : Real
81 0 : NearestReporterCoordinatesFunction::timeDerivative(Real t, const Point & p) const
82 : {
83 0 : const std::array<std::pair<Real, std::size_t>, 2> tv = findNearestPoint(t, p);
84 :
85 0 : if (tv[0].second == tv[1].second)
86 : return 0.0;
87 :
88 0 : const Real told = tv[0].first;
89 0 : const Real tnew = tv[1].first;
90 0 : const Real vold = _values[tv[0].second];
91 0 : const Real vnew = _values[tv[1].second];
92 0 : return (vnew - vold) / (tnew - told);
93 : }
94 :
95 : std::vector<Real>
96 2223168 : NearestReporterCoordinatesFunction::parameterGradient(Real t, const Point & p) const
97 : {
98 2223168 : const std::array<std::pair<Real, std::size_t>, 2> tv = findNearestPoint(t, p);
99 2223168 : std::vector<Real> param_grad(_nval, 0.0);
100 :
101 2223168 : if (tv[0].second == tv[1].second)
102 2219712 : param_grad[tv[0].second] = 1;
103 : else
104 : {
105 3456 : const Real told = tv[0].first;
106 3456 : const Real tnew = tv[1].first;
107 3456 : param_grad[tv[0].second] = (tnew - t) / (tnew - told);
108 3456 : param_grad[tv[1].second] = (t - told) / (tnew - told);
109 : }
110 2223168 : return param_grad;
111 : }
112 :
113 : void
114 238 : NearestReporterCoordinatesFunction::buildCoordinateMapping() const
115 : {
116 : // Do some size checks
117 238 : _nval = std::max({_coordx.size(), _coordy.size(), _coordz.size(), _coordt.size()});
118 238 : if (_nval == 0)
119 0 : paramError("value", "At least one coordinate vector must not be empty.");
120 238 : else if (!_coordx.empty() && _coordx.size() != _nval)
121 3 : paramError("x_coord_name",
122 : "Number of x coordinates (",
123 : _coordx.size(),
124 : ") does not match number of values (",
125 : _nval,
126 : ").");
127 235 : else if (!_coordy.empty() && _coordy.size() != _nval)
128 3 : paramError("y_coord_name",
129 : "Number of y coordinates (",
130 : _coordy.size(),
131 : ") does not match number of values (",
132 : _nval,
133 : ").");
134 232 : else if (!_coordz.empty() && _coordz.size() != _nval)
135 3 : paramError("z_coord_name",
136 : "Number of z coordinates (",
137 : _coordz.size(),
138 : ") does not match number of values (",
139 : _nval,
140 : ").");
141 229 : else if (!_coordt.empty() && _coordt.size() != _nval)
142 3 : paramError("time_name",
143 : "Number of times (",
144 : _coordt.size(),
145 : ") does not match number of values (",
146 : _nval,
147 : ").");
148 :
149 : // Find times for each unique coordinate
150 : _coord_mapping.clear();
151 1752 : for (const auto & i : make_range(_nval))
152 : {
153 : Point pt;
154 1526 : pt(0) = _coordx.empty() ? 0.0 : _coordx[i];
155 1526 : pt(1) = _coordy.empty() ? 0.0 : _coordy[i];
156 1526 : pt(2) = _coordz.empty() ? 0.0 : _coordz[i];
157 1526 : const Real time = _coordt.empty() ? 0.0 : _coordt[i];
158 :
159 : std::vector<std::pair<Real, std::size_t>> * vec = nullptr;
160 5557 : for (auto & it : _coord_mapping)
161 4641 : if (pt.absolute_fuzzy_equals(it.first))
162 : {
163 610 : vec = &it.second;
164 : break;
165 : }
166 : if (!vec)
167 916 : vec = &_coord_mapping[pt];
168 :
169 1526 : vec->emplace_back(time, i);
170 1526 : std::sort(vec->begin(),
171 : vec->end(),
172 : [](const std::pair<Real, Real> & a, const std::pair<Real, Real> & b)
173 : { return a.first < b.first; });
174 : }
175 226 : }
176 :
177 : std::array<std::pair<Real, std::size_t>, 2>
178 48335247 : NearestReporterCoordinatesFunction::findNearestPoint(Real t, const Point & p) const
179 : {
180 48335247 : if (_coord_mapping.empty())
181 238 : buildCoordinateMapping();
182 :
183 : // Make sure values is correct size
184 48335235 : if (_values.size() != _nval)
185 3 : paramError("value_name",
186 : "Size of value vector (",
187 : _values.size(),
188 : ") does not match number of coordinates specified (",
189 : _nval,
190 : ").");
191 :
192 : const auto & tval =
193 48335232 : std::min_element(_coord_mapping.begin(),
194 : _coord_mapping.end(),
195 137006784 : [&p](const std::pair<Point, std::vector<std::pair<Real, std::size_t>>> & p1,
196 : const std::pair<Point, std::vector<std::pair<Real, std::size_t>>> & p2)
197 137006784 : { return (p - p1.first).norm_sq() < (p - p2.first).norm_sq(); })
198 : ->second;
199 :
200 48335232 : if (tval.size() == 1 || MooseUtils::absoluteFuzzyLessEqual(t, tval[0].first))
201 48325632 : return {tval[0], tval[0]};
202 9600 : else if (MooseUtils::absoluteFuzzyGreaterEqual(t, tval.back().first))
203 960 : return {tval.back(), tval.back()};
204 : else
205 12480 : for (std::size_t ti = 1; ti < tval.size(); ++ti)
206 12480 : if (MooseUtils::absoluteFuzzyGreaterEqual(tval[ti].first, t))
207 8640 : return {tval[ti - 1], tval[ti]};
208 :
209 0 : mooseError("Internal error: unable to find nearest point.");
210 : return std::array<std::pair<Real, std::size_t>, 2>();
211 : }
|
