Skip to content

Commit 3629862

Browse files
crswong888crswong888
committed
Create a function to perform baseline correction on acceleration time histories (refs idaholab#296)
1 parent 502e5e1 commit 3629862

File tree

4 files changed

+559
-0
lines changed

4 files changed

+559
-0
lines changed

‎include/functions/BaselineCorrection.h

+75
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
/*************************************************/
2+
/* DO NOT MODIFY THIS HEADER */
3+
/* */
4+
/* MASTODON */
5+
/* */
6+
/* (c) 2015 Battelle Energy Alliance, LLC */
7+
/* ALL RIGHTS RESERVED */
8+
/* */
9+
/* Prepared by Battelle Energy Alliance, LLC */
10+
/* With the U. S. Department of Energy */
11+
/* */
12+
/* See COPYRIGHT for full restrictions */
13+
/*************************************************/
14+
15+
// This code was implemented in collaboration with Christopher J. Wong
16+
// ([email protected]) from the University of Utah.
17+
18+
#pragma once
19+
20+
// MOOSE includes
21+
#include "Function.h"
22+
#include "LinearInterpolation.h"
23+
24+
// Forward Declarations
25+
class BaselineCorrection;
26+
27+
/**
28+
* Applies a baseline correction to an accceleration time history using least
29+
* squares polynomial fits and outputs the adjusted acceleration
30+
*/
31+
class BaselineCorrection : public Function
32+
{
33+
public:
34+
static InputParameters validParams();
35+
36+
BaselineCorrection(const InputParameters & parameters);
37+
38+
virtual Real value(Real t, const Point & /*P*/) const override;
39+
40+
protected:
41+
/// Newmark integration parameters
42+
const Real & _gamma;
43+
const Real & _beta;
44+
45+
/// set which kinematic variables a polynomial fit will be applied to
46+
const bool _fit_accel;
47+
const bool _fit_vel;
48+
const bool _fit_disp;
49+
50+
/// order used for the least squares polynomial fit
51+
const unsigned int _order;
52+
53+
/// acceleration time history variables from input
54+
std::vector<Real> _time;
55+
std::vector<Real> _accel;
56+
57+
/// adjusted (corrected) acceleration ordinates
58+
std::vector<Real> _adj_accel;
59+
60+
/// object to output linearly interpolated corrected acceleration ordinates
61+
std::unique_ptr<LinearInterpolation> _linear_interp;
62+
63+
/// function value scale factor
64+
const Real & _scale_factor;
65+
66+
private:
67+
/// Applies baseline correction to raw acceleration time history
68+
void applyCorrection();
69+
70+
/// Reads and builds data from supplied CSV file
71+
void buildFromFile();
72+
73+
/// Builds data from pairs of `time_values` and `acceleration_values'
74+
void buildFromXandY();
75+
};

‎include/utils/BaselineCorrectionUtils.h

+72
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
/*************************************************/
2+
/* DO NOT MODIFY THIS HEADER */
3+
/* */
4+
/* MASTODON */
5+
/* */
6+
/* (c) 2015 Battelle Energy Alliance, LLC */
7+
/* ALL RIGHTS RESERVED */
8+
/* */
9+
/* Prepared by Battelle Energy Alliance, LLC */
10+
/* With the U. S. Department of Energy */
11+
/* */
12+
/* See COPYRIGHT for full restrictions */
13+
/*************************************************/
14+
15+
// This code was implemented in collaboration with Christopher J. Wong
16+
// ([email protected]) from the University of Utah.
17+
18+
#pragma once
19+
20+
// LIBMESH includes
21+
#include "DenseMatrix.h"
22+
#include "libmesh/dense_vector.h"
23+
24+
/**
25+
* This namespace contains the functions used for the calculations corresponding
26+
* to the time history adjustment procedure in BaselineCorrection
27+
**/
28+
namespace BaselineCorrectionUtils
29+
{
30+
/// Evaluates an integral over a single time step with Newmark-beta method
31+
/// Also is used as simple trapezoidal rule when gamma = 0.5.
32+
Real newmarkGammaIntegrate(const Real & u_ddot_old,
33+
const Real & u_ddot,
34+
const Real & u_dot_old,
35+
const Real & gamma,
36+
const Real & dt);
37+
38+
/// Evaluates a double integral over a single time step with Newmark-beta method
39+
Real newmarkBetaIntegrate(const Real & u_ddot_old,
40+
const Real & u_ddot,
41+
const Real & u_dot_old,
42+
const Real & u_old,
43+
const Real & beta,
44+
const Real & dt);
45+
46+
/// Solves linear normal equation for minimum acceleration square error
47+
DenseVector<Real> getAccelerationFitCoeffs(unsigned int order,
48+
const std::vector<Real> & accel,
49+
const std::vector<Real> & t,
50+
const unsigned int & num_steps,
51+
const Real & gamma);
52+
53+
/// Solves linear normal equation for minimum velocity square error
54+
DenseVector<Real> getVelocityFitCoeffs(unsigned int order,
55+
const std::vector<Real> & accel,
56+
const std::vector<Real> & vel,
57+
const std::vector<Real> & t,
58+
const unsigned int & num_steps,
59+
const Real & beta);
60+
61+
/// Solves linear normal equation for minimum displacement square error
62+
DenseVector<Real> getDisplacementFitCoeffs(unsigned int order,
63+
const std::vector<Real> & disp,
64+
const std::vector<Real> & t,
65+
const unsigned int & num_steps);
66+
67+
/// Evaluates the least squares polynomials over at a single time step
68+
std::vector<Real> computePolynomials(unsigned int order,
69+
const DenseVector<Real> & coeffs,
70+
const Real & t);
71+
72+
} // namespace BaselineCorrectionUtils

‎src/functions/BaselineCorrection.C

+234
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,234 @@
1+
/*************************************************/
2+
/* DO NOT MODIFY THIS HEADER */
3+
/* */
4+
/* MASTODON */
5+
/* */
6+
/* (c) 2015 Battelle Energy Alliance, LLC */
7+
/* ALL RIGHTS RESERVED */
8+
/* */
9+
/* Prepared by Battelle Energy Alliance, LLC */
10+
/* With the U. S. Department of Energy */
11+
/* */
12+
/* See COPYRIGHT for full restrictions */
13+
/*************************************************/
14+
15+
// This code was implemented in collaboration with Christopher J. Wong
16+
// ([email protected]) from the University of Utah.
17+
18+
// MASTODON includes
19+
#include "BaselineCorrection.h"
20+
#include "BaselineCorrectionUtils.h"
21+
22+
// MOOSE includes
23+
#include "DelimitedFileReader.h"
24+
25+
registerMooseObject("MastodonApp", BaselineCorrection);
26+
27+
InputParameters
28+
BaselineCorrection::validParams()
29+
{
30+
InputParameters params = Function::validParams();
31+
32+
params.addParam<FileName>("data_file",
33+
"The name of a CSV file containing raw acceleration time history data.");
34+
params.addParam<std::string>("time_name",
35+
"The header name of the column which contains the time values in the data file. If not "
36+
"specified, they are assumed to be in the first column index.");
37+
params.addParam<std::string>("acceleration_name",
38+
"The header name for the column which contains the acceleration values in the data file. If "
39+
"not specified, they are assumed to be in the second column index.");
40+
params.addParam<std::vector<Real>>("time_values", "The time abscissa values.");
41+
params.addParam<std::vector<Real>>("acceleration_values", "The acceleration ordinate values.");
42+
43+
params.addRequiredParam<Real>("gamma", "The gamma parameter for Newmark time integration.");
44+
params.addRequiredParam<Real>("beta", "The beta parameter for Newmark time integration.");
45+
46+
params.addParam<bool>("fit_acceleration", true,
47+
"If set to \"true\", the acceleration time history will be adjusted using a polynomial fit "
48+
"of the acceleration data.");
49+
params.addParam<bool>("fit_velocity", false,
50+
"If set to \"true\", the acceleration time history will be adjusted using a polynomial fit "
51+
"of the velocity data obtained by integration.");
52+
params.addParam<bool>("fit_displacement", false,
53+
"If set to \"true\", the acceleration time history will be adjusted using a polynomial fit "
54+
"of the displacement data obtained by double-integration.");
55+
params.addRequiredRangeCheckedParam<unsigned int>("order", "(0 < order) & (order < 10)",
56+
"The order of the polynomial fit(s) used to adjust the nominal time histories (coefficients "
57+
"of higher order polynomials can be difficult to compute and the method generally becomes "
58+
"unstable when order >= 10).");
59+
60+
params.addParam<Real>("scale_factor", 1.0,
61+
"A scale factor to be applied to the adjusted acceleration time history.");
62+
params.declareControllable("scale_factor");
63+
64+
return params;
65+
}
66+
67+
BaselineCorrection::BaselineCorrection(const InputParameters & parameters)
68+
: Function(parameters),
69+
_gamma(getParam<Real>("gamma")),
70+
_beta(getParam<Real>("beta")),
71+
_fit_accel(getParam<bool>("fit_acceleration")),
72+
_fit_vel(getParam<bool>("fit_velocity")),
73+
_fit_disp(getParam<bool>("fit_displacement")),
74+
_order(getParam<unsigned int>("order")),
75+
_scale_factor(getParam<Real>("scale_factor"))
76+
{
77+
// determine data source and check parameter consistency
78+
if (isParamValid("data_file") && !isParamValid("time_values") &&
79+
!isParamValid("acceleration_values"))
80+
buildFromFile();
81+
else if (!isParamValid("data_file") && isParamValid("time_values") &&
82+
isParamValid("acceleration_values"))
83+
buildFromXandY();
84+
else
85+
mooseError("In BaselineCorrection ",
86+
_name,
87+
": Either `data_file` or `time_values` and `acceleration_values` must be specified "
88+
"exclusively.");
89+
90+
// size checks
91+
if (_time.size() != _accel.size())
92+
mooseError("In BaselineCorrection ",
93+
_name,
94+
": The length of time and acceleration data must be equal.");
95+
if (_time.size() == 0 || _accel.size() == 0)
96+
mooseError("In BaselineCorrection ",
97+
_name,
98+
": The length of time and acceleration data must be > 0.");
99+
100+
// ensure that at least one best-fit will be created
101+
if (!_fit_accel && !_fit_vel && !_fit_disp)
102+
mooseWarning("Warning in " + name() +
103+
". Computation of a polynomial fit is set to \"false\" for all three "
104+
"kinematic variables. No adjustments will occur and the output will be the "
105+
"raw acceleration time history.");
106+
107+
// apply baseline correction to raw acceleration time history
108+
applyCorrection();
109+
110+
// try building a linear interpolation object
111+
try
112+
{
113+
_linear_interp = libmesh_make_unique<LinearInterpolation>(_time, _adj_accel);
114+
}
115+
catch (std::domain_error & e)
116+
{
117+
mooseError("In BaselineCorrection ", _name, ": ", e.what());
118+
}
119+
}
120+
121+
Real
122+
BaselineCorrection::value(Real t, const Point & /*P*/) const
123+
{
124+
return _scale_factor * _linear_interp->sample(t);
125+
}
126+
127+
void
128+
BaselineCorrection::applyCorrection()
129+
{
130+
// store a reference to final array index
131+
unsigned int index_end = _accel.size() - 1;
132+
133+
// Compute unadjusted velocity and displacment time histories
134+
Real dt;
135+
std::vector<Real> vel, disp;
136+
vel.push_back(0); disp.push_back(0);
137+
for (unsigned int i = 0; i < index_end; ++i)
138+
{
139+
dt = _time[i+1] - _time[i];
140+
141+
vel.push_back(BaselineCorrectionUtils::newmarkGammaIntegrate(
142+
_accel[i], _accel[i+1], vel[i], _gamma, dt));
143+
disp.push_back(BaselineCorrectionUtils::newmarkBetaIntegrate(
144+
_accel[i], _accel[i+1], vel[i], disp[i], _beta, dt));
145+
}
146+
147+
// initialize polyfits and adjusted time history arrays as the nominal ones
148+
DenseVector<Real> coeffs;
149+
_adj_accel = _accel;
150+
std::vector<Real> p_fit, adj_vel = vel, adj_disp = disp;
151+
152+
// adjust time histories with acceleration fit if desired
153+
if (_fit_accel)
154+
{
155+
coeffs = BaselineCorrectionUtils::getAccelerationFitCoeffs(
156+
_order, _adj_accel, _time, index_end, _gamma);
157+
158+
for (unsigned int i = 0; i <= index_end; ++i)
159+
{
160+
p_fit = BaselineCorrectionUtils::computePolynomials(_order, coeffs, _time[i]);
161+
_adj_accel[i] -= p_fit[0];
162+
adj_vel[i] -= p_fit[1];
163+
adj_disp[i] -= p_fit[2];
164+
}
165+
}
166+
167+
// adjust with velocity fit
168+
if (_fit_vel)
169+
{
170+
coeffs = BaselineCorrectionUtils::getVelocityFitCoeffs(
171+
_order, _adj_accel, adj_vel, _time, index_end, _beta);
172+
173+
for (unsigned int i = 0; i <= index_end; ++i)
174+
{
175+
p_fit = BaselineCorrectionUtils::computePolynomials(_order, coeffs, _time[i]);
176+
_adj_accel[i] -= p_fit[0];
177+
adj_disp[i] -= p_fit[2];
178+
}
179+
}
180+
181+
// adjust with displacement fit
182+
if (_fit_disp)
183+
{
184+
coeffs = BaselineCorrectionUtils::getDisplacementFitCoeffs(_order, adj_disp, _time, index_end);
185+
186+
for (unsigned int i = 0; i <= index_end; ++i)
187+
{
188+
p_fit = BaselineCorrectionUtils::computePolynomials(_order, coeffs, _time[i]);
189+
_adj_accel[i] -= p_fit[0];
190+
}
191+
}
192+
}
193+
194+
void
195+
BaselineCorrection::buildFromFile()
196+
{
197+
// Read data from CSV file
198+
MooseUtils::DelimitedFileReader reader(getParam<FileName>("data_file"), &_communicator);
199+
reader.read();
200+
201+
// Check if specific column headers were input
202+
const bool time_header = isParamValid("time_name"),
203+
accel_header = isParamValid("acceleration_name");
204+
205+
if (time_header && !accel_header)
206+
mooseError("In BaselineCorrection ",
207+
_name,
208+
": A column header name was specified for the for the time data. Please specify a ",
209+
"header for the acceleration data using the 'accelertation_name' parameter.");
210+
else if (!time_header && accel_header)
211+
mooseError("In BaselineCorrection ",
212+
_name,
213+
": A column header name was specified for the for the acceleration data. Please "
214+
"specify a header for the time data using the 'time_name' parameter.");
215+
216+
// Obtain acceleration time history from file data
217+
if (time_header && accel_header)
218+
{
219+
_time = reader.getData(getParam<std::string>("time_name"));
220+
_accel = reader.getData(getParam<std::string>("acceleration_name"));
221+
}
222+
else
223+
{
224+
_time = reader.getData(0);
225+
_accel = reader.getData(1);
226+
}
227+
}
228+
229+
void
230+
BaselineCorrection::buildFromXandY()
231+
{
232+
_time = getParam<std::vector<Real>>("time_values");
233+
_accel = getParam<std::vector<Real>>("acceleration_values");
234+
}

0 commit comments

Comments
 (0)