diff --git a/include/functions/BaselineCorrection.h b/include/functions/BaselineCorrection.h new file mode 100644 index 0000000000..0ed11fe6b6 --- /dev/null +++ b/include/functions/BaselineCorrection.h @@ -0,0 +1,75 @@ +/*************************************************/ +/* 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 + */ +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; + + /// set which kinematic variables a polynomial fit will be applied to + const bool _fit_accel; + const bool _fit_vel; + const bool _fit_disp; + + /// order used for the least squares polynomial fit + const 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..724869e511 --- /dev/null +++ b/include/utils/BaselineCorrectionUtils.h @@ -0,0 +1,72 @@ +/*************************************************/ +/* 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..426df01bd0 --- /dev/null +++ b/src/functions/BaselineCorrection.C @@ -0,0 +1,234 @@ +/*************************************************/ +/* 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.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.addParam("fit_acceleration", true, + "If set to \"true\", the acceleration time history will be adjusted using a polynomial fit " + "of the acceleration data."); + params.addParam("fit_velocity", false, + "If set to \"true\", the acceleration time history will be adjusted using a polynomial fit " + "of the velocity data obtained by integration."); + params.addParam("fit_displacement", false, + "If set to \"true\", the acceleration time history will be adjusted using a polynomial fit " + "of the displacement data obtained by double-integration."); + params.addRequiredRangeCheckedParam("order", "(0 < order) & (order < 10)", + "The order of the polynomial fit(s) used to adjust the nominal time histories (coefficients " + "of higher order polynomials can be difficult to compute and the method generally becomes " + "unstable when order >= 10)."); + + 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")), + _fit_accel(getParam("fit_acceleration")), + _fit_vel(getParam("fit_velocity")), + _fit_disp(getParam("fit_displacement")), + _order(getParam("order")), + _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."); + + // ensure that at least one best-fit will be created + if (!_fit_accel && !_fit_vel && !_fit_disp) + mooseWarning("Warning in " + name() + + ". Computation of a polynomial fit is set to \"false\" for all three " + "kinematic variables. No adjustments will occur and the output will be the " + "raw acceleration time history."); + + // 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 (_fit_accel) + { + 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 (_fit_vel) + { + 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 (_fit_disp) + { + 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 'accelertation_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..02d9a5fe08 --- /dev/null +++ b/src/utils/BaselineCorrectionUtils.C @@ -0,0 +1,178 @@ +/*************************************************/ +/* 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; +}