[WIP] Implementing time-history baseline correction (closes idaholab#296

)

doc/content/source/vectorpostprocessors/LeastSquaresBaselineCorrection.md

@@ -0,0 +1,26 @@
+# LeastSquaresBaselineCorrection (DRAFT)
+
+The `LeastSquaresBaselineCorrection` applies a baseline correction to an acceleration time-history, $A(t)$, contained in some other VectorPostprocessor object, e.g., `ResponseHistoryBuilder` or `CSVReader`. The corrections are performed by first computing the nominal velocity, $V(t)$, and displacement, $D(t)$, time-histories by the [Newmark-beta method](manuals/theory/index.md#time-integration). An $n$-th order polynomial approximation of the original, or *uncorrected*, acceleration, velocity, or displacement is then found by the method of least squares and is subtracted from the uncorrected acceleration to obtain the corrected acceleration, $\tilde{A}$. That is, the corrected acceleration time-history is given by the following:
+
+\begin{equation}
+\label{corrected-accel}
+\tilde{A}(t)=A(t)-P_{A}(t)-P_{V}(t)-P_{D}(t)
+\end{equation}
+
+Here, $P_{A}$, $P_{V}$, or $P_{D}$ are the least squares approximations of the uncorrected acceleration, velocity, and displacement time histories, respectively. [corrected-accel] is then integrated to obtain the corrected velocity, $\tilde{V}$, and displacement, $\tilde{D}$. This object, ultimately, stores the corrected time-histories as VectorPostprocessor data. Note that it is possible to apply corrections using a best-fit of only one or any combination of the three kinematic variables (see [#usage]), for which, any unused polynomials would result in a zero value in [corrected-accel].
+
+## Usage id=usage
+
+The object requires, as input, some sort of record of an acceleration time-history stored in a MOOSE VectorPostprocessor object specified with the `vectorpostprocessor` parameter. The VectorPostprocessor columns which store the abscissa (time) values and the ordinate (acceleration) values are then specified with `time_name` and `acceleration_name` parameters, respectively. This data could be obtained from any number of sources, for example, it could be the nodal acceleration data obtained from the finite element solution itself, or it could be the raw data from an accelelogram. If the case is the latter example, this object could be useful for using the corrected time-history data elsewhere in the input file. Note if the corrections are to be used as inputs for the simulation, the user must take care to set the `execute_on`, `force_preaux`, or `force_preic` parameters accordingly.
+
+The order of the approximating polynomials, $n$, is specified with the `order` parameter and can take on a maximum integer value of 9. The corrections are highly sensitive to the order, and so it is recommended that a user run with several trial values until the desired result is achieved. The `beta` and `gamma` parameters are those used for Newmark Integration, and, ideally, should correspond to the values used elsewhere throughout the simulation. Finally, the `fit_acceleration`, `fit_velocity`, and `fit_displacement` variables control whether or not to use the polynomials described in [corrected-accel] when formulating the corrected time-histories.
+
+## Example Input Syntax
+
+!listing test/tests/vectorpostprocessors/least_squares_baseline_correction/least_squares_baseline_correction.i block=VectorPostprocessors
+
+!syntax parameters /VectorPostprocessors/LeastSquaresBaselineCorrection
+
+!syntax inputs /VectorPostprocessors/LeastSquaresBaselineCorrection
+
+!syntax children /VectorPostprocessors/LeastSquaresBaselineCorrection

include/utils/BaselineCorrectionUtils.h

@@ -0,0 +1,66 @@
+// This code was implemented in collaboration with Christopher J. Wong
+// ([email protected]) from the University of Utah.
+
+/**
+ * This namespace contains the functions used for the calculations corresponding
+ * to the time history adjustment procedure in LeastSquaresBaselineCorrection
+ **/
+
+#ifndef BASELINECORRECTIONUTILS_H
+#define BASELINECORRECTIONUTILS_H
+
+// MOOSE includes
+#include "DenseMatrix.h"
+
+// LIBMESH includes
+#include "libmesh/dense_vector.h"
+
+// Forward declarations
+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<Real> getAccelerationFitCoeffs(unsigned int order,
+                                           const std::vector<Real> & accel,
+                                           const std::vector<Real> & t,
+                                           const unsigned int & num_steps,
+                                           const Real & gamma);
+
+// Solves linear normal equation for minimum velocity square error
+DenseVector<Real> getVelocityFitCoeffs(unsigned int order,
+                                       const std::vector<Real> & accel,
+                                       const std::vector<Real> & vel,
+                                       const std::vector<Real> & t,
+                                       const unsigned int & num_steps,
+                                       const Real & beta);
+
+// Solves linear normal equation for minimum displacement square error
+DenseVector<Real> getDisplacementFitCoeffs(unsigned int order,
+                                           const std::vector<Real> & disp,
+                                           const std::vector<Real> & t,
+                                           const unsigned int & num_steps);
+
+// Evaluates the least squares polynomials over at a single time step
+std::vector<Real> computePolynomials(unsigned int order,
+                                     const DenseVector<Real> & coeffs,
+                                     const Real & t);
+
+}
+
+#endif // BASELINECORRECTIONUTILS_H

include/vectorpostprocessors/LeastSquaresBaselineCorrection.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       */
+/*************************************************/
+
+#ifndef LEASTSQUARESBASELINECORRECTION_H
+#define LEASTSQUARESBASELINECORRECTION_H
+
+// MOOSE includes
+#include "GeneralVectorPostprocessor.h"
+
+// Forward Declarations
+class LeastSquaresBaselineCorrection;
+
+template <>
+InputParameters validParams<LeastSquaresBaselineCorrection>();
+
+/**
+ * Applies a baseline correction to an accceleration time history contained in another
+ * vectorpostprocessor using least squares polynomial fits and outputs the adjusted
+ * acceleration, velocity, and displacement time histories
+ */
+class LeastSquaresBaselineCorrection : public GeneralVectorPostprocessor
+{
+public:
+  LeastSquaresBaselineCorrection(const InputParameters & parameters);
+  virtual void initialize() override;
+  virtual void execute() override;
+
+protected:
+  // acceleration time history variables from specified vectorpostprocessor
+  const VectorPostprocessorValue & _accel;
+  const VectorPostprocessorValue & _t;
+
+  // order used for the least squares polynomial fit
+  const unsigned int _order;
+
+  // 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;
+
+  // the variables used to write out the adjusted time histories
+  VectorPostprocessorValue & _time;
+  VectorPostprocessorValue & _adj_accel;
+  VectorPostprocessorValue & _adj_vel;
+  VectorPostprocessorValue & _adj_disp;
+
+  // wether to output nominal time histories from Newmark integration alone
+  const bool _out_unadj;
+
+  // the variables used to write out the nominal time histories (if requested)
+  VectorPostprocessorValue * _unadj_accel;
+  VectorPostprocessorValue * _unadj_vel;
+  VectorPostprocessorValue * _unadj_disp;
+};
+
+#endif // LEASTSQUARESBASELINECORRECTION_H

src/utils/BaselineCorrectionUtils.C

@@ -0,0 +1,163 @@
+// This code was implemented in collaboration with Christopher J. Wong
+// ([email protected]) 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<Real>
+BaselineCorrectionUtils::getAccelerationFitCoeffs(unsigned int order,
+                                                  const std::vector<Real> & accel,
+                                                  const std::vector<Real> & t,
+                                                  const unsigned int & num_steps,
+                                                  const Real & gamma)
+{
+  unsigned int num_rows = order + 1; /* no. of eqns to solve for coefficients */
+  DenseMatrix<Real> mat(num_rows, num_rows);
+  DenseVector<Real> rhs(num_rows);
+  DenseVector<Real> 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<Real>
+BaselineCorrectionUtils::getVelocityFitCoeffs(unsigned int order,
+                                              const std::vector<Real> & accel,
+                                              const std::vector<Real> & vel,
+                                              const std::vector<Real> & t,
+                                              const unsigned int & num_steps,
+                                              const Real & beta)
+{
+  unsigned int num_rows = order + 1; /* no. of eqns to solve for coefficients */
+  DenseMatrix<Real> mat(num_rows, num_rows);
+  DenseVector<Real> rhs(num_rows);
+  DenseVector<Real> 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<Real>
+BaselineCorrectionUtils::getDisplacementFitCoeffs(unsigned int order,
+                                                  const std::vector<Real> & disp,
+                                                  const std::vector<Real> & t,
+                                                  const unsigned int & num_steps)
+{
+  unsigned int num_rows = order + 1;
+  DenseMatrix<Real> mat(num_rows, num_rows);
+  DenseVector<Real> rhs(num_rows);
+  DenseVector<Real> 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<Real>
+BaselineCorrectionUtils::computePolynomials(unsigned int order,
+                                            const DenseVector<Real> & coeffs,
+                                            const Real & t)
+{
+  std::vector<Real> p_fit(3); /* accel polyfit and its derivatives */
+  for (unsigned int k = 0; k < order + 1; ++k) /* compute polynomials */
+  {
+    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;
+}