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 180 : GeneralSensorPostprocessor::validParams()
24 : {
25 180 : InputParameters params = GeneralPostprocessor::validParams();
26 360 : params.addRequiredParam<PostprocessorName>("input_signal", "The input signal postprocessor");
27 360 : params.addParam<FunctionName>(
28 360 : "drift_function", 0.0, "Drift function describing signal drift over time");
29 360 : params.addParam<Real>("vector_size", 1000000000.0, "The maximum size of vector to be saved");
30 180 : params.addClassDescription("This is a GeneralSensorPostprocessor, and functions as a base class "
31 : "for other sensor postprocessors");
32 360 : params.addParam<FunctionName>(
33 360 : "efficiency_function", 0.8, "Efficiency function describing efficiency over time");
34 360 : params.addParam<FunctionName>(
35 : "noise_std_dev_function",
36 360 : 0.5,
37 : "Standard deviation of noise function describing noise standard deviation over time");
38 360 : params.addParam<FunctionName>(
39 360 : "signalToNoise_function", 0.2, "Signal to noise ratio for the modeled sensor");
40 360 : params.addParam<FunctionName>(
41 360 : "delay_function", 0.1, "Delay function describing the delay over time");
42 360 : params.addParam<FunctionName>(
43 : "uncertainty_std_dev_function",
44 360 : 0.1,
45 : "Standard deviation of uncertainty function describing uncertainty std dev over time");
46 360 : params.addParam<FunctionName>("R_function", 1, "The function R for the integration term.");
47 360 : params.addParam<Real>("proportional_weight", 1, "The weight assigned to the proportional term");
48 360 : params.addParam<Real>("integral_weight", 0, "The weight assigned to the integral term");
49 360 : params.addParam<unsigned int>("seed", 1, "Seed for the random number generator");
50 180 : return params;
51 0 : }
52 :
53 90 : GeneralSensorPostprocessor::GeneralSensorPostprocessor(const InputParameters & parameters)
54 : : GeneralPostprocessor(parameters),
55 90 : _input_signal(getPostprocessorValue("input_signal")),
56 90 : _drift_function(getFunction("drift_function")),
57 180 : _vector_size(getParam<Real>("vector_size")),
58 90 : _pp_old(getPostprocessorValueOld("input_signal")),
59 90 : _efficiency_function(getFunction("efficiency_function")),
60 90 : _noise_std_dev_function(getFunction("noise_std_dev_function")),
61 90 : _delay_function(getFunction("delay_function")),
62 90 : _signalToNoise_function(getFunction("signalToNoise_function")),
63 90 : _uncertainty_std_dev_function(getFunction("uncertainty_std_dev_function")),
64 90 : _R_function(getFunction("R_function")),
65 180 : _proportional_weight(getParam<Real>("proportional_weight")),
66 180 : _integral_weight(getParam<Real>("integral_weight")),
67 180 : _time_values(declareRestartableData<std::vector<Real>>("time_values")),
68 180 : _input_signal_values(declareRestartableData<std::vector<Real>>("input_signal_values")),
69 180 : _integrand(declareRestartableData<std::vector<Real>>("integrand")),
70 180 : _R_function_values(declareRestartableData<std::vector<Real>>("R_function_values")),
71 180 : _t_step_old(declareRestartableData<int>("gsp_t_step_old", -1)),
72 180 : _seed(getParam<unsigned int>("seed")),
73 90 : _delay_value(0) // Initialize delay_value at 0
74 : {
75 90 : }
76 :
77 : void
78 3928 : GeneralSensorPostprocessor::initialize()
79 : {
80 3928 : Real drift_value = _drift_function.value(_t);
81 3928 : Real efficiency_value = _efficiency_function.value(_t);
82 3928 : Real signalToNoise_value = _signalToNoise_function.value(_t);
83 : // setting seed for random number generator
84 3928 : _rng.seed(_seed);
85 3928 : Real noise_std_dev = _noise_std_dev_function.value(_t);
86 : Real noise_value = _rng.randNormal(0, noise_std_dev);
87 3928 : 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 3928 : _delay_value = _delay_function.value(_t);
91 :
92 : // if the problem is steady-state
93 3928 : if (!_fe_problem.isTransient())
94 : {
95 : // output
96 13 : _sensor_value = drift_value +
97 13 : 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 3915 : 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 3915 : _time_values.push_back(_t);
116 3915 : _input_signal_values.push_back(_input_signal);
117 :
118 : // Check if the size is greater than _vector_size
119 3915 : 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 3915 : Real _input_signal_delayed = getDelayedInputSignal();
126 :
127 : // computing integral term
128 3915 : Real term_for_integration = _input_signal + signalToNoise_value * noise_value;
129 3915 : _integrand.push_back(term_for_integration);
130 3915 : _integration_value = getIntegral(_integrand);
131 :
132 : // output
133 3915 : Real proportional_value = _input_signal_delayed + signalToNoise_value * noise_value;
134 3915 : _sensor_value = drift_value +
135 3915 : efficiency_value * (_proportional_weight * proportional_value +
136 3915 : _integral_weight * _integration_value) +
137 : uncertainty_value;
138 :
139 : // Update old time step
140 3915 : _t_step_old = _t_step;
141 : }
142 3928 : }
143 :
144 : PostprocessorValue
145 3928 : GeneralSensorPostprocessor::getValue() const
146 : {
147 3928 : return _sensor_value;
148 : }
149 :
150 : std::vector<Real>
151 4566 : 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 4566 : result.reserve(vec1.size());
160 350880 : for (const auto i : index_range(vec1))
161 346314 : result.push_back(vec1[i] * vec2[i]);
162 4566 : return result;
163 0 : }
164 :
165 : Real
166 4566 : GeneralSensorPostprocessor::getDelayedInputSignal()
167 : {
168 : // delayed input signal
169 4566 : Real t_desired = _t - _delay_value;
170 : Real input_signal_delayed;
171 :
172 4566 : if (t_desired < _time_values[0])
173 : input_signal_delayed = 0;
174 :
175 : // Modified cases where desired time isin the time value list
176 4540 : else if (t_desired == _time_values[0])
177 26 : input_signal_delayed = _input_signal_values[0];
178 :
179 : // linear interpolation
180 4514 : else if (t_desired > _time_values[0] && t_desired <= _t && t_desired >= _t - _dt)
181 2550 : input_signal_delayed = _input_signal + (t_desired - _t) * (_pp_old - _input_signal) / (-_dt);
182 :
183 1964 : else if (t_desired > _time_values[0] && t_desired < _t - _dt)
184 : {
185 1964 : LinearInterpolation time_and_input_signal(_time_values, _input_signal_values);
186 1964 : input_signal_delayed = time_and_input_signal.sample(t_desired);
187 1964 : }
188 :
189 : else
190 0 : mooseError("Unhandled case in interpolation values.");
191 :
192 4566 : return input_signal_delayed;
193 : }
194 :
195 : vector<Real>
196 3915 : GeneralSensorPostprocessor::getRVector()
197 : {
198 3915 : _R_function_values.push_back(_R_function.value(_t));
199 3915 : return _R_function_values;
200 : }
201 :
202 : Real
203 4566 : GeneralSensorPostprocessor::getIntegral(std::vector<Real> integrand)
204 : {
205 : Real integration_value;
206 4566 : getRVector();
207 : // calculate the total integrand vector including previous integrand and exponential terms
208 9132 : 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 4566 : if (_time_values.size() > 2)
213 : {
214 4488 : 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 4488 : if (_time_values.size() < 30)
220 : n = 30;
221 : else
222 : {
223 3434 : n = _time_values.size();
224 : // if interval is not a multiple of 3, make it
225 6855 : while (static_cast<int>(n) % 3 != 0)
226 3421 : n = n + 1;
227 : }
228 :
229 4488 : 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 373341 : for (Real i = 0; i < n + 1; i++)
236 : {
237 368853 : Real new_time = _time_values[0] + i * h;
238 368853 : time_vec_simp.push_back(new_time);
239 368853 : Real new_integrand = spline.SplineInterpolation::sample(new_time);
240 368853 : integrand_vec_simp.push_back(new_integrand);
241 : }
242 :
243 : // initialize integral with endpoints
244 4488 : integration_value = integrand_vec_simp[0] + integrand_vec_simp[n];
245 :
246 : // Sum for points not at endpoints
247 364365 : for (int i = 1; i < n; ++i)
248 : {
249 359877 : if (i % 3 == 0)
250 : // Every third point
251 116967 : integration_value += 2 * integrand_vec_simp[i];
252 : else
253 242910 : integration_value += 3 * integrand_vec_simp[i];
254 : }
255 :
256 : // Multiply by 3h/8
257 4488 : integration_value *= (3 * h) / 8.0;
258 4488 : }
259 :
260 : else
261 : {
262 : // performing integration by trapezoidal rule, not accurate
263 : // LinearInterpolation object with two vectors
264 78 : LinearInterpolation integral(_time_values, integrand);
265 : // integrate the vectors
266 78 : integration_value = integral.integrate();
267 78 : }
268 :
269 4566 : return integration_value;
270 : }
|
