From c86182e8ebc782e8f1d4f9b7e47aaabb48423eab Mon Sep 17 00:00:00 2001 From: crswong888 Date: Fri, 10 Jul 2020 14:39:36 -0600 Subject: [PATCH] Create a function to perform baseline correction on acceleration time histories (closes #296) --- include/functions/BaselineCorrection.h | 70 +++++++ include/utils/BaselineCorrectionUtils.h | 71 +++++++ src/functions/BaselineCorrection.C | 247 ++++++++++++++++++++++++ src/utils/BaselineCorrectionUtils.C | 180 +++++++++++++++++ 4 files changed, 568 insertions(+) create mode 100644 include/functions/BaselineCorrection.h create mode 100644 include/utils/BaselineCorrectionUtils.h create mode 100644 src/functions/BaselineCorrection.C create mode 100644 src/utils/BaselineCorrectionUtils.C diff --git a/include/functions/BaselineCorrection.h b/include/functions/BaselineCorrection.h new file mode 100644 index 0000000000..56933d26e7 --- /dev/null +++ b/include/functions/BaselineCorrection.h @@ -0,0 +1,70 @@ +/*************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MASTODON */ +/* */ +/* (c) 2015 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/*************************************************/ + +// This code was implemented in collaboration with Christopher J. Wong +// (chris.wong@utah.edu) from the University of Utah. + +#pragma once + +// MOOSE includes +#include "Function.h" +#include "LinearInterpolation.h" + +// Forward Declarations +class BaselineCorrection; + +/** + * Applies a baseline correction to an accceleration time history using least squares polynomial + * fits and outputs the adjusted acceleration with linear interpolation. + */ +class BaselineCorrection : public Function +{ +public: + static InputParameters validParams(); + + BaselineCorrection(const InputParameters & parameters); + + virtual Real value(Real t, const Point & /*P*/) const override; + +protected: + /// Newmark integration parameters + const Real & _gamma; + const Real & _beta; + + /// order used for the least squares polynomial fit + unsigned int _order; + + /// acceleration time history variables from input + std::vector _time; + std::vector _accel; + + /// adjusted (corrected) acceleration ordinates + std::vector _adj_accel; + + /// object to output linearly interpolated corrected acceleration ordinates + std::unique_ptr _linear_interp; + + /// function value scale factor + const Real & _scale_factor; + +private: + /// Applies baseline correction to raw acceleration time history + void applyCorrection(); + + /// Reads and builds data from supplied CSV file + void buildFromFile(); + + /// Builds data from pairs of `time_values` and `acceleration_values' + void buildFromXandY(); +}; diff --git a/include/utils/BaselineCorrectionUtils.h b/include/utils/BaselineCorrectionUtils.h new file mode 100644 index 0000000000..fcc02f5d56 --- /dev/null +++ b/include/utils/BaselineCorrectionUtils.h @@ -0,0 +1,71 @@ +/*************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MASTODON */ +/* */ +/* (c) 2015 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/*************************************************/ + +// This code was implemented in collaboration with Christopher J. Wong +// (chris.wong@utah.edu) from the University of Utah. + +#pragma once + +// LIBMESH includes +#include "DenseMatrix.h" +#include "libmesh/dense_vector.h" + +/** + * This namespace contains the functions used for the calculations corresponding + * to the time history adjustment procedure in BaselineCorrection + **/ +namespace BaselineCorrectionUtils +{ +/// Evaluates an integral over a single time step with Newmark-beta method +/// Also is used as simple trapezoidal rule when gamma = 0.5. +Real newmarkGammaIntegrate(const Real & u_ddot_old, + const Real & u_ddot, + const Real & u_dot_old, + const Real & gamma, + const Real & dt); + +/// Evaluates a double integral over a single time step with Newmark-beta method +Real newmarkBetaIntegrate(const Real & u_ddot_old, + const Real & u_ddot, + const Real & u_dot_old, + const Real & u_old, + const Real & beta, + const Real & dt); + +/// Solves linear normal equation for minimum acceleration square error +DenseVector getAccelerationFitCoeffs(unsigned int order, + const std::vector & accel, + const std::vector & t, + const unsigned int & num_steps, + const Real & gamma); + +/// Solves linear normal equation for minimum velocity square error +DenseVector getVelocityFitCoeffs(unsigned int order, + const std::vector & accel, + const std::vector & vel, + const std::vector & t, + const unsigned int & num_steps, + const Real & beta); + +/// Solves linear normal equation for minimum displacement square error +DenseVector getDisplacementFitCoeffs(unsigned int order, + const std::vector & disp, + const std::vector & t, + const unsigned int & num_steps); + +/// Evaluates the least squares polynomials over at a single time step +std::vector +computePolynomials(unsigned int order, const DenseVector & coeffs, const Real & t); + +} // namespace BaselineCorrectionUtils diff --git a/src/functions/BaselineCorrection.C b/src/functions/BaselineCorrection.C new file mode 100644 index 0000000000..43d3adf09d --- /dev/null +++ b/src/functions/BaselineCorrection.C @@ -0,0 +1,247 @@ +/*************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MASTODON */ +/* */ +/* (c) 2015 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/*************************************************/ + +// This code was implemented in collaboration with Christopher J. Wong +// (chris.wong@utah.edu) from the University of Utah. + +// MASTODON includes +#include "BaselineCorrection.h" +#include "BaselineCorrectionUtils.h" + +// MOOSE includes +#include "DelimitedFileReader.h" + +registerMooseObject("MastodonApp", BaselineCorrection); + +InputParameters +BaselineCorrection::validParams() +{ + InputParameters params = Function::validParams(); + params.addClassDescription("Applies a baseline correction to an accceleration time history " + "using least squares polynomial fits and outputs the adjusted " + "acceleration with linear interpolation."); + + params.addParam( + "data_file", "The name of a CSV file containing raw acceleration time history data."); + params.addParam( + "time_name", + "The header name of the column which contains the time values in the data file. If not " + "specified, they are assumed to be in the first column index."); + params.addParam( + "acceleration_name", + "The header name for the column which contains the acceleration values in the data file. If " + "not specified, they are assumed to be in the second column index."); + + params.addParam>("time_values", "The time abscissa values."); + params.addParam>("acceleration_values", "The acceleration ordinate values."); + + params.addRequiredParam("gamma", "The gamma parameter for Newmark time integration."); + params.addRequiredParam("beta", "The beta parameter for Newmark time integration."); + + params.addRangeCheckedParam( + "accel_fit_order", + "(0 <= accel_fit_order) & (accel_fit_order < 10)", + "If this is provided, the acceleration time history will be adjusted using an nth-order " + "polynomial fit of the nominal acceleration data, where n = accel_fit_order (only integer " + "values from 0 to 9 are supported)."); + params.addRangeCheckedParam( + "vel_fit_order", + "(0 <= vel_fit_order) & (vel_fit_order < 10)", + "If this is provided, the acceleration time history will be adjusted using an nth-order " + "polynomial fit of the nominal velocity data, where n = vel_fit_order (only integer values " + "from 0 to 9 are supported)."); + params.addRangeCheckedParam( + "disp_fit_order", + "(0 <= disp_fit_order) & (disp_fit_order < 10)", + "If this is provided, the acceleration time history will be adjusted using an nth-order " + "polynomial fit of the nominal displacement data, where n = disp_fit_order (only integer " + "values from 0 to 9 are supported)."); + + params.addParam("scale_factor", + 1.0, + "A scale factor to be applied to the adjusted acceleration time history."); + params.declareControllable("scale_factor"); + + return params; +} + +BaselineCorrection::BaselineCorrection(const InputParameters & parameters) + : Function(parameters), + _gamma(getParam("gamma")), + _beta(getParam("beta")), + _scale_factor(getParam("scale_factor")) +{ + // determine data source and check parameter consistency + if (isParamValid("data_file") && !isParamValid("time_values") && + !isParamValid("acceleration_values")) + buildFromFile(); + else if (!isParamValid("data_file") && isParamValid("time_values") && + isParamValid("acceleration_values")) + buildFromXandY(); + else + mooseError("In BaselineCorrection ", + _name, + ": Either `data_file` or `time_values` and `acceleration_values` must be specified " + "exclusively."); + + // size checks + if (_time.size() != _accel.size()) + mooseError("In BaselineCorrection ", + _name, + ": The length of time and acceleration data must be equal."); + if (_time.size() == 0 || _accel.size() == 0) + mooseError( + "In BaselineCorrection ", _name, ": The length of time and acceleration data must be > 0."); + + // check that at lease one least squares fit will be applied + if (!isParamValid("accel_fit_order") && !isParamValid("vel_fit_order") && + !isParamValid("disp_fit_order")) + mooseError("In BaselineCorrection ", + _name, + ": No values were input for parameters 'accel_fit_order', 'vel_fit_order', or " + "'disp_fit_order'. Please specify an integer value from 0 to 9 for at least one " + "of these parameters."); + + // apply baseline correction to raw acceleration time history + applyCorrection(); + + // try building a linear interpolation object + try + { + _linear_interp = libmesh_make_unique(_time, _adj_accel); + } + catch (std::domain_error & e) + { + mooseError("In BaselineCorrection ", _name, ": ", e.what()); + } +} + +Real +BaselineCorrection::value(Real t, const Point & /*P*/) const +{ + return _scale_factor * _linear_interp->sample(t); +} + +void +BaselineCorrection::applyCorrection() +{ + // store a reference to final array index + unsigned int index_end = _accel.size() - 1; + + // Compute unadjusted velocity and displacment time histories + Real dt; + std::vector vel, disp; + vel.push_back(0); + disp.push_back(0); + for (unsigned int i = 0; i < index_end; ++i) + { + dt = _time[i + 1] - _time[i]; + + vel.push_back(BaselineCorrectionUtils::newmarkGammaIntegrate( + _accel[i], _accel[i + 1], vel[i], _gamma, dt)); + disp.push_back(BaselineCorrectionUtils::newmarkBetaIntegrate( + _accel[i], _accel[i + 1], vel[i], disp[i], _beta, dt)); + } + + // initialize polyfits and adjusted time history arrays as the nominal ones + DenseVector coeffs; + _adj_accel = _accel; + std::vector p_fit, adj_vel = vel, adj_disp = disp; + + // adjust time histories with acceleration fit if desired + if (isParamValid("accel_fit_order")) + { + _order = getParam("accel_fit_order"); + coeffs = BaselineCorrectionUtils::getAccelerationFitCoeffs( + _order, _adj_accel, _time, index_end, _gamma); + + for (unsigned int i = 0; i <= index_end; ++i) + { + p_fit = BaselineCorrectionUtils::computePolynomials(_order, coeffs, _time[i]); + _adj_accel[i] -= p_fit[0]; + adj_vel[i] -= p_fit[1]; + adj_disp[i] -= p_fit[2]; + } + } + + // adjust with velocity fit + if (isParamValid("vel_fit_order")) + { + _order = getParam("vel_fit_order"); + coeffs = BaselineCorrectionUtils::getVelocityFitCoeffs( + _order, _adj_accel, adj_vel, _time, index_end, _beta); + + for (unsigned int i = 0; i <= index_end; ++i) + { + p_fit = BaselineCorrectionUtils::computePolynomials(_order, coeffs, _time[i]); + _adj_accel[i] -= p_fit[0]; + adj_disp[i] -= p_fit[2]; + } + } + + // adjust with displacement fit + if (isParamValid("disp_fit_order")) + { + _order = getParam("disp_fit_order"); + coeffs = BaselineCorrectionUtils::getDisplacementFitCoeffs(_order, adj_disp, _time, index_end); + + for (unsigned int i = 0; i <= index_end; ++i) + { + p_fit = BaselineCorrectionUtils::computePolynomials(_order, coeffs, _time[i]); + _adj_accel[i] -= p_fit[0]; + } + } +} + +void +BaselineCorrection::buildFromFile() +{ + // Read data from CSV file + MooseUtils::DelimitedFileReader reader(getParam("data_file"), &_communicator); + reader.read(); + + // Check if specific column headers were input + const bool time_header = isParamValid("time_name"), + accel_header = isParamValid("acceleration_name"); + + if (time_header && !accel_header) + mooseError("In BaselineCorrection ", + _name, + ": A column header name was specified for the for the time data. Please specify a ", + "header for the acceleration data using the 'acceleration_name' parameter."); + else if (!time_header && accel_header) + mooseError("In BaselineCorrection ", + _name, + ": A column header name was specified for the for the acceleration data. Please " + "specify a header for the time data using the 'time_name' parameter."); + + // Obtain acceleration time history from file data + if (time_header && accel_header) + { + _time = reader.getData(getParam("time_name")); + _accel = reader.getData(getParam("acceleration_name")); + } + else + { + _time = reader.getData(0); + _accel = reader.getData(1); + } +} + +void +BaselineCorrection::buildFromXandY() +{ + _time = getParam>("time_values"); + _accel = getParam>("acceleration_values"); +} diff --git a/src/utils/BaselineCorrectionUtils.C b/src/utils/BaselineCorrectionUtils.C new file mode 100644 index 0000000000..148788f57f --- /dev/null +++ b/src/utils/BaselineCorrectionUtils.C @@ -0,0 +1,180 @@ +/*************************************************/ +/* DO NOT MODIFY THIS HEADER */ +/* */ +/* MASTODON */ +/* */ +/* (c) 2015 Battelle Energy Alliance, LLC */ +/* ALL RIGHTS RESERVED */ +/* */ +/* Prepared by Battelle Energy Alliance, LLC */ +/* With the U. S. Department of Energy */ +/* */ +/* See COPYRIGHT for full restrictions */ +/*************************************************/ + +// This code was implemented in collaboration with Christopher J. Wong +// (chris.wong@utah.edu) from the University of Utah. + +// MASTODON includes +#include "BaselineCorrectionUtils.h" + +Real +BaselineCorrectionUtils::newmarkGammaIntegrate(const Real & u_ddot_old, + const Real & u_ddot, + const Real & u_dot_old, + const Real & gamma, + const Real & dt) +{ + return u_dot_old + (1 - gamma) * dt * u_ddot_old + gamma * dt * u_ddot; +} + +Real +BaselineCorrectionUtils::newmarkBetaIntegrate(const Real & u_ddot_old, + const Real & u_ddot, + const Real & u_dot_old, + const Real & u_old, + const Real & beta, + const Real & dt) +{ + return u_old + dt * u_dot_old + (0.5 - beta) * dt * dt * u_ddot_old + beta * dt * dt * u_ddot; +} + +DenseVector +BaselineCorrectionUtils::getAccelerationFitCoeffs(unsigned int order, + const std::vector & accel, + const std::vector & t, + const unsigned int & num_steps, + const Real & gamma) +{ + unsigned int num_rows = order + 1; + DenseMatrix mat(num_rows, num_rows); + DenseVector rhs(num_rows); + DenseVector coeffs(num_rows); + + // compute matrix of linear normal equation + for (unsigned int row = 0; row < num_rows; ++row) + { + for (unsigned int col = 0; col < num_rows; ++col) + { + mat(row, col) = + pow(t[t.size() - 1], row + col + 1) * (col * col + 3 * col + 2) / (row + col + 1); + } + } + + // compute vector of integrals on right-hand side of linear normal equation + Real dt, u_ddot_old, u_ddot; + for (unsigned int i = 0; i < num_steps; ++i) + { + dt = t[i + 1] - t[i]; + for (unsigned int row = 0; row < num_rows; ++row) + { + u_ddot_old = pow(t[i], row) * accel[i]; + u_ddot = pow(t[i + 1], row) * accel[i + 1]; + + rhs(row) += newmarkGammaIntegrate(u_ddot_old, u_ddot, 0.0, gamma, dt); + } + } + + // solve the system using libMesh lu factorization + mat.lu_solve(rhs, coeffs); + return coeffs; +} + +DenseVector +BaselineCorrectionUtils::getVelocityFitCoeffs(unsigned int order, + const std::vector & accel, + const std::vector & vel, + const std::vector & t, + const unsigned int & num_steps, + const Real & beta) +{ + unsigned int num_rows = order + 1; + DenseMatrix mat(num_rows, num_rows); + DenseVector rhs(num_rows); + DenseVector coeffs(num_rows); + + // compute matrix of linear normal equation + for (unsigned int row = 0; row < num_rows; ++row) + { + for (unsigned int col = 0; col < num_rows; ++col) + { + mat(row, col) = pow(t[t.size() - 1], row + col + 3) * (col + 2) / (row + col + 3); + } + } + + // compute vector of integrals on right-hand side of linear normal equation + Real dt, u_ddot_old, u_ddot, u_dot_old; + for (unsigned int i = 0; i < num_steps; ++i) + { + dt = t[i + 1] - t[i]; + for (unsigned int row = 0; row < num_rows; ++row) + { + u_dot_old = pow(t[i], row + 1) * vel[i]; + u_ddot_old = pow(t[i], row + 1) * accel[i] + (row + 1) * pow(t[i], row) * vel[i]; + u_ddot = pow(t[i + 1], row + 1) * accel[i + 1] + (row + 1) * pow(t[i + 1], row) * vel[i + 1]; + + rhs(row) += newmarkBetaIntegrate(u_ddot_old, u_ddot, u_dot_old, 0.0, beta, dt); + } + } + + // solve the system using libMesh lu factorization + mat.lu_solve(rhs, coeffs); + return coeffs; +} + +DenseVector +BaselineCorrectionUtils::getDisplacementFitCoeffs(unsigned int order, + const std::vector & disp, + const std::vector & t, + const unsigned int & num_steps) +{ + unsigned int num_rows = order + 1; + DenseMatrix mat(num_rows, num_rows); + DenseVector rhs(num_rows); + DenseVector coeffs(num_rows); + + // computer matrix of linear normal equation + for (unsigned int row = 0; row < num_rows; ++row) + { + for (unsigned int col = 0; col < num_rows; ++col) + { + mat(row, col) = pow(t[t.size() - 1], row + col + 5) / (row + col + 5); + } + } + + // compute vector of integrals on right-hand side of linear normal equation + Real dt, u_old, u; + for (unsigned int i = 0; i < num_steps; ++i) + { + dt = t[i + 1] - t[i]; + for (unsigned int row = 0; row < num_rows; ++row) + { + u_old = pow(t[i], row + 2) * disp[i]; + u = pow(t[i + 1], row + 2) * disp[i + 1]; + + // note: newmarkGamma with gamma = 0.5 is trapezoidal rule + rhs(row) += newmarkGammaIntegrate(u_old, u, 0.0, 0.5, dt); + } + } + + // solve the system using libMesh lu factorization + mat.lu_solve(rhs, coeffs); + return coeffs; +} + +std::vector +BaselineCorrectionUtils::computePolynomials(unsigned int order, + const DenseVector & coeffs, + const Real & t) +{ + // Compute the best-fit polynomial of the acceleration and its derivatives + std::vector p_fit(3); + for (unsigned int k = 0; k < order + 1; ++k) + { + p_fit[0] += (k * k + 3 * k + 2) * coeffs(k) * pow(t, k); + p_fit[1] += (k + 2) * coeffs(k) * pow(t, k + 1); + p_fit[2] += coeffs(k) * pow(t, k + 2); + } + + return p_fit; +}