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 : #include <vector>
10 : #include <cmath>
11 : #include <algorithm>
12 : #include <numeric>
13 :
14 : #include "GeneralSensorPostprocessor.h"
15 : #include "Function.h"
16 : #include "MooseRandom.h"
17 : #include "LinearInterpolation.h"
18 : #include "SplineInterpolation.h"
19 :
20 : registerMooseObject("MiscApp", GeneralSensorPostprocessor);
21 :
22 : InputParameters
23 80 : GeneralSensorPostprocessor::validParams()
24 : {
25 80 : InputParameters params = GeneralPostprocessor::validParams();
26 160 : params.addRequiredParam<PostprocessorName>("input_signal", "The input signal postprocessor");
27 160 : params.addParam<FunctionName>(
28 160 : "drift_function", 0.0, "Drift function describing signal drift over time");
29 160 : params.addParam<Real>("vector_size", 1000000000.0, "The maximum size of vector to be saved");
30 80 : params.addClassDescription("This is a GeneralSensorPostprocessor, and functions as a base class "
31 : "for other sensor postprocessors");
32 160 : params.addParam<FunctionName>(
33 160 : "efficiency_function", 0.8, "Efficiency function describing efficiency over time");
34 160 : params.addParam<FunctionName>(
35 : "noise_std_dev_function",
36 160 : 0.5,
37 : "Standard deviation of noise function describing noise standard deviation over time");
38 160 : params.addParam<FunctionName>(
39 160 : "signalToNoise_function", 0.2, "Signal to noise ratio for the modeled sensor");
40 160 : params.addParam<FunctionName>(
41 160 : "delay_function", 0.1, "Delay function describing the delay over time");
42 160 : params.addParam<FunctionName>(
43 : "uncertainty_std_dev_function",
44 160 : 0.1,
45 : "Standard deviation of uncertainty function describing uncertainty std dev over time");
46 160 : params.addParam<FunctionName>("R_function", 1, "The function R for the integration term.");
47 160 : params.addParam<Real>("proportional_weight", 1, "The weight assigned to the proportional term");
48 160 : params.addParam<Real>("integral_weight", 0, "The weight assigned to the integral term");
49 160 : params.addParam<unsigned int>("seed", 1, "Seed for the random number generator");
50 80 : return params;
51 0 : }
52 :
53 40 : GeneralSensorPostprocessor::GeneralSensorPostprocessor(const InputParameters & parameters)
54 : : GeneralPostprocessor(parameters),
55 40 : _input_signal(getPostprocessorValue("input_signal")),
56 40 : _drift_function(getFunction("drift_function")),
57 80 : _vector_size(getParam<Real>("vector_size")),
58 40 : _pp_old(getPostprocessorValueOld("input_signal")),
59 40 : _efficiency_function(getFunction("efficiency_function")),
60 40 : _noise_std_dev_function(getFunction("noise_std_dev_function")),
61 40 : _delay_function(getFunction("delay_function")),
62 40 : _signalToNoise_function(getFunction("signalToNoise_function")),
63 40 : _uncertainty_std_dev_function(getFunction("uncertainty_std_dev_function")),
64 40 : _R_function(getFunction("R_function")),
65 80 : _proportional_weight(getParam<Real>("proportional_weight")),
66 80 : _integral_weight(getParam<Real>("integral_weight")),
67 80 : _time_values(declareRestartableData<std::vector<Real>>("time_values")),
68 80 : _input_signal_values(declareRestartableData<std::vector<Real>>("input_signal_values")),
69 80 : _integrand(declareRestartableData<std::vector<Real>>("integrand")),
70 80 : _R_function_values(declareRestartableData<std::vector<Real>>("R_function_values")),
71 80 : _t_step_old(declareRestartableData<int>("gsp_t_step_old", -1)),
72 80 : _seed(getParam<unsigned int>("seed")),
73 40 : _delay_value(0) // Initialize delay_value at 0
74 : {
75 40 : }
76 :
77 : void
78 2116 : GeneralSensorPostprocessor::initialize()
79 : {
80 2116 : Real drift_value = _drift_function.value(_t);
81 2116 : Real efficiency_value = _efficiency_function.value(_t);
82 2116 : Real signalToNoise_value = _signalToNoise_function.value(_t);
83 : // setting seed for random number generator
84 2116 : _rng.seed(_seed);
85 2116 : Real noise_std_dev = _noise_std_dev_function.value(_t);
86 : Real noise_value = _rng.randNormal(0, noise_std_dev);
87 2116 : Real uncertainty_std_dev = _uncertainty_std_dev_function.value(_t);
88 : Real uncertainty_value = _rng.randNormal(0, uncertainty_std_dev);
89 : // Added line to calculate delay_value from the delay function
90 2116 : _delay_value = _delay_function.value(_t);
91 :
92 : // if the problem is steady-state
93 2116 : if (!_fe_problem.isTransient())
94 : {
95 : // output
96 7 : _sensor_value = drift_value +
97 7 : efficiency_value * (_input_signal + signalToNoise_value * noise_value) +
98 : uncertainty_value;
99 : }
100 :
101 : // if the problem is transient
102 : else
103 : {
104 : // Remove last element if we are repeating the timestep
105 : mooseAssert(_t_step_old <= _t_step,
106 : "The old time step needs to be behind or the same as the current time step.");
107 2109 : if (_t_step_old == _t_step)
108 : {
109 2 : _time_values.pop_back();
110 2 : _input_signal_values.pop_back();
111 2 : _integrand.pop_back();
112 2 : _R_function_values.pop_back();
113 : }
114 :
115 2109 : _time_values.push_back(_t);
116 2109 : _input_signal_values.push_back(_input_signal);
117 :
118 : // Check if the size is greater than _vector_size
119 2109 : if (_time_values.size() > _vector_size && _input_signal_values.size() > _vector_size)
120 : {
121 : // Remove the first 10 elements if size is larger
122 : _time_values.erase(_time_values.begin(), _time_values.begin() + 10);
123 0 : _input_signal_values.erase(_input_signal_values.begin(), _input_signal_values.begin() + 10);
124 : }
125 2109 : Real _input_signal_delayed = getDelayedInputSignal();
126 :
127 : // computing integral term
128 2109 : Real term_for_integration = _input_signal + signalToNoise_value * noise_value;
129 2109 : _integrand.push_back(term_for_integration);
130 2109 : _integration_value = getIntegral(_integrand);
131 :
132 : // output
133 2109 : Real proportional_value = _input_signal_delayed + signalToNoise_value * noise_value;
134 2109 : _sensor_value = drift_value +
135 2109 : efficiency_value * (_proportional_weight * proportional_value +
136 2109 : _integral_weight * _integration_value) +
137 : uncertainty_value;
138 :
139 : // Update old time step
140 2109 : _t_step_old = _t_step;
141 : }
142 2116 : }
143 :
144 : PostprocessorValue
145 2116 : GeneralSensorPostprocessor::getValue() const
146 : {
147 2116 : return _sensor_value;
148 : }
149 :
150 : std::vector<Real>
151 2460 : GeneralSensorPostprocessor::elementwiseMultiply(std::vector<Real> & vec1, std::vector<Real> & vec2)
152 : {
153 : // Ensure both vectors have the same size
154 : mooseAssert(vec1.size() == vec2.size(),
155 : "Vectors must have the same size for element-wise multiplication.");
156 : // Create a result vector
157 : std::vector<Real> result;
158 : // Perform element-wise multiplication
159 2460 : result.reserve(vec1.size());
160 189018 : for (const auto i : index_range(vec1))
161 186558 : result.push_back(vec1[i] * vec2[i]);
162 2460 : return result;
163 0 : }
164 :
165 : Real
166 2460 : GeneralSensorPostprocessor::getDelayedInputSignal()
167 : {
168 : // delayed input signal
169 2460 : Real t_desired = _t - _delay_value;
170 : Real input_signal_delayed;
171 :
172 2460 : if (t_desired < _time_values[0])
173 : input_signal_delayed = 0;
174 :
175 : // Modified cases where desired time isin the time value list
176 2446 : else if (t_desired == _time_values[0])
177 14 : input_signal_delayed = _input_signal_values[0];
178 :
179 : // linear interpolation
180 2432 : else if (t_desired > _time_values[0] && t_desired <= _t && t_desired >= _t - _dt)
181 1374 : input_signal_delayed = _input_signal + (t_desired - _t) * (_pp_old - _input_signal) / (-_dt);
182 :
183 1058 : else if (t_desired > _time_values[0] && t_desired < _t - _dt)
184 : {
185 1058 : LinearInterpolation time_and_input_signal(_time_values, _input_signal_values);
186 1058 : input_signal_delayed = time_and_input_signal.sample(t_desired);
187 1058 : }
188 :
189 : else
190 0 : mooseError("Unhandled case in interpolation values.");
191 :
192 2460 : return input_signal_delayed;
193 : }
194 :
195 : vector<Real>
196 2109 : GeneralSensorPostprocessor::getRVector()
197 : {
198 2109 : _R_function_values.push_back(_R_function.value(_t));
199 2109 : return _R_function_values;
200 : }
201 :
202 : Real
203 2460 : GeneralSensorPostprocessor::getIntegral(std::vector<Real> integrand)
204 : {
205 : Real integration_value;
206 2460 : getRVector();
207 : // calculate the total integrand vector including previous integrand and exponential terms
208 4920 : integrand = elementwiseMultiply(integrand, _R_function_values);
209 :
210 : // if there are more than two datapoints, do simpson's 3/8 rule for integration
211 : // else do trapezoidal rule
212 2460 : if (_time_values.size() > 2)
213 : {
214 2418 : SplineInterpolation spline(_time_values, integrand);
215 : // number of intervals
216 : Real n;
217 :
218 : // if number of datapoints is less than 30, use spline interpolation to get 30 intervals
219 2418 : if (_time_values.size() < 30)
220 : n = 30;
221 : else
222 : {
223 1850 : n = _time_values.size();
224 : // if interval is not a multiple of 3, make it
225 3693 : while (static_cast<int>(n) % 3 != 0)
226 1843 : n = n + 1;
227 : }
228 :
229 2418 : Real h = (_time_values.back() - _time_values[0]) / n; // distance between time values
230 : // time vector for simpson integration
231 : vector<Real> time_vec_simp;
232 : // integrand vector for simpson integration
233 : vector<Real> integrand_vec_simp;
234 :
235 201117 : for (Real i = 0; i < n + 1; i++)
236 : {
237 198699 : Real new_time = _time_values[0] + i * h;
238 198699 : time_vec_simp.push_back(new_time);
239 198699 : Real new_integrand = spline.SplineInterpolation::sample(new_time);
240 198699 : integrand_vec_simp.push_back(new_integrand);
241 : }
242 :
243 : // initialize integral with endpoints
244 2418 : integration_value = integrand_vec_simp[0] + integrand_vec_simp[n];
245 :
246 : // Sum for points not at endpoints
247 196281 : for (int i = 1; i < n; ++i)
248 : {
249 193863 : if (i % 3 == 0)
250 : // Every third point
251 63009 : integration_value += 2 * integrand_vec_simp[i];
252 : else
253 130854 : integration_value += 3 * integrand_vec_simp[i];
254 : }
255 :
256 : // Multiply by 3h/8
257 2418 : integration_value *= (3 * h) / 8.0;
258 2418 : }
259 :
260 : else
261 : {
262 : // performing integration by trapezoidal rule, not accurate
263 : // LinearInterpolation object with two vectors
264 42 : LinearInterpolation integral(_time_values, integrand);
265 : // integrate the vectors
266 42 : integration_value = integral.integrate();
267 42 : }
268 :
269 2460 : return integration_value;
270 : }
|
