Merge pull request idaholab#269 from cbolisetti/rsp

Adding period column to the response spectra calculator VPP. Closes #…

doc/content/source/vectorpostprocessors/ResponseSpectraCalculator.md

@@ -1,21 +1,21 @@
 # ResponseSpectraCalculator
 
-The `ResponseSpectraCalculator` calculates the acceleration, velocity and displacement response spectra for a set of response histories calculated by the `ResponseHistoryBuilder` VectorPostprocessor. The response spectra will be calculated as functions of frequency and the outputs will be stored in a csv file.
+The `ResponseSpectraCalculator` calculates the acceleration, velocity and displacement response spectra for a set of response histories calculated by the `ResponseHistoryBuilder` VectorPostprocessor. The response spectra will be calculated as functions of frequency and period, and the outputs will be stored in a csv file.
 
 ## Usage
 
 For example, the input file below contains a `ResponseHistoryBuilder` VectorPostprocessor called 'accel_hist', which builds the response histories of the acceleration variables, accel_x and accel_y at nodes 8 and 9. It is desired that the corresponding response spectra are calculated. This can be calculated by `ResponseSpectraCalculator` using the 'accel_spec' block as follows.
 
 !listing test/tests/vectorpostprocessors/response_spectra_calculator/response_spectra_calculator.i block=VectorPostprocessors
 
-!alert note
-`ResponseSpectraCalculator` assumes that the response histories calculated in the vectorpostprocessor, provided in the `vectorpostprocessor` parameter in the input are acceleration histories, and will not perform any checks to ensure that they are accelerations. The user should ensure this by defining the corresponding `ResponseHistoryBuilder` appropriately. 
+!alert note title=Acceleration histories only
+`ResponseSpectraCalculator` assumes that the response histories provided through the `vectorpostprocessor` input parameter (e.g., in the 'accel_hist' block above) are acceleration histories. No checks are performed to ensure that they are accelerations. The user should ensure this by defining the corresponding `ResponseHistoryBuilder` block appropriately.
 
 Special care must be taken to ensure that the response spectra calculations are performed only at the end of the simulation (last time step) and not at every time step. In order to accomplish this, an 'Outputs' block must be created as shown below. This block (named 'out') ensures that the CSV output is only calculated at the final timestep.
 
 !listing test/tests/vectorpostprocessors/response_spectra_calculator/response_spectra_calculator.i block=Outputs
 
-The name of the Outputs block 'out' should be provided in the response spectra VectorPostprocessor block as `outputs = out` as shown above. For this input, a csv file of the response spectra is created, as shown below. Note that the columns in the csv file are organized according to the node numbers and variable names. For example, the spectral accelerations for accel_x at node 8 are under node_8_accel_x_sd, node_8_accel_x_sv, and node_8_accel_x_sa, which are the displacement, velocity and acceleration spectra, respectively.
+The name of the Outputs block 'out' should be provided in the response spectra VectorPostprocessor block as `outputs = out` as shown above. For this input, a csv file of the response spectra is created, as shown below. Note that the columns in the csv file are organized according to the node numbers and variable names. For example, the spectral accelerations for accel_x at node 8 are under node_8_accel_x_sd, node_8_accel_x_sv, and node_8_accel_x_sa, which are the displacement, velocity and acceleration spectra, respectively. Note that the corresponding frequency and period columns are also included in the output file.
 
 !listing test/tests/vectorpostprocessors/response_spectra_calculator/gold/response_spectra_calculator_out_accel_spec.csv
 

include/vectorpostprocessors/ResponseSpectraCalculator.h

@@ -61,6 +61,9 @@ class ResponseSpectraCalculator : public GeneralVectorPostprocessor
   /// Reference to the frequency vector.
   VectorPostprocessorValue & _frequency;
 
+  /// Reference to the period vector.
+  VectorPostprocessorValue & _period;
+
   /// Vector of pointers to the vectors of spectral values.
   std::vector<VectorPostprocessorValue *> _spectrum;
 

src/utils/MastodonUtils.C

@@ -22,14 +22,16 @@ MastodonUtils::responseSpectrum(const Real & freq_start,
                                 const Real & xi,
                                 const Real & reg_dt)
 {
-  std::vector<Real> freq_vec, aspec_vec, vspec_vec, dspec_vec;
+  std::vector<Real> freq_vec, per_vec, aspec_vec, vspec_vec, dspec_vec;
   Real logdf, om_n, om_d, dt2, dis1, vel1, acc1, dis2, vel2, acc2, pdmax, kd;
   for (std::size_t n = 0; n < freq_num; ++n)
   {
-    // Building the frequency vector. Frequencies are distributed
-    // uniformly in the log scale.
+    // Building the frequency vector and the period vector.
+    // Frequencies are distributed uniformly in the log scale.
+    // Periods are calculated as inverse of frequency
     logdf = (std::log10(freq_end) - std::log10(freq_start)) / (freq_num - 1);
     freq_vec.push_back(pow(10.0, std::log10(freq_start) + n * logdf));
+    per_vec.push_back(1.0 / freq_vec[n]);
     om_n = 2.0 * 3.141593 * freq_vec[n]; // om_n = 2*pi*f
     om_d = om_n * xi;
     dis1 = 0.0;
@@ -55,7 +57,7 @@ MastodonUtils::responseSpectrum(const Real & freq_start,
     vspec_vec.push_back(pdmax * om_n);
     aspec_vec.push_back(pdmax * om_n * om_n);
   }
-  return {freq_vec, dspec_vec, vspec_vec, aspec_vec};
+  return {freq_vec, per_vec, dspec_vec, vspec_vec, aspec_vec};
 }
 
 std::vector<std::vector<Real>>

src/vectorpostprocessors/Fragility.C

@@ -237,7 +237,7 @@ Fragility::calcDemandsFromFile(unsigned int bin)
           demand_sample, demand_time, _dtsim)[1]; // regularize the demand sample
       demand_sample_spectrum =
           MastodonUtils::responseSpectrum(0.01, 100, 401, demand_sample, _ssc_xi, _dtsim);
-      LinearInterpolation spectraldemand(demand_sample_spectrum[0], demand_sample_spectrum[3]);
+      LinearInterpolation spectraldemand(demand_sample_spectrum[0], demand_sample_spectrum[4]);
       stoc_demands[k] = spectraldemand.sample(_ssc_freq);
       k++;
     }

src/vectorpostprocessors/HousnerSpectrumIntensity.C

@@ -84,7 +84,6 @@ HousnerSpectrumIntensity::initialize()
 void
 HousnerSpectrumIntensity::execute()
 {
-  std::vector<Real> frequency;
   std::vector<Real> period;
   std::vector<Real> vel_spectrum;
   Real freq_start = 1 / _per_end;
@@ -99,11 +98,8 @@ HousnerSpectrumIntensity::execute()
     // Calculation of the response spectrum.
     std::vector<std::vector<Real>> var_spectrum = MastodonUtils::responseSpectrum(
         freq_start, freq_end, _per_num, reg_vector[1], _xi, _reg_dt);
-    frequency = var_spectrum[0];
-    vel_spectrum = var_spectrum[2];
-    period.resize(frequency.size());
-    for (unsigned int j = 0; j < frequency.size(); ++j)
-      period[j] = 1 / frequency[j];
+    period = var_spectrum[1];
+    vel_spectrum = var_spectrum[3];
     std::reverse(period.begin(), period.end());
     std::reverse(vel_spectrum.begin(), vel_spectrum.end());
     LinearInterpolation hsi_calc(period, vel_spectrum);

src/vectorpostprocessors/ResponseSpectraCalculator.C

@@ -47,6 +47,7 @@ ResponseSpectraCalculator::ResponseSpectraCalculator(const InputParameters & par
     _freq_num(getParam<unsigned int>("num_frequencies")),
     _reg_dt(getParam<Real>("regularize_dt")),
     _frequency(declareVector("frequency")),
+    _period(declareVector("period")),
     // Time vector from the response history builder vector postprocessor
     _history_time(getVectorPostprocessorValue("vectorpostprocessor", "time"))
 
@@ -55,7 +56,9 @@ ResponseSpectraCalculator::ResponseSpectraCalculator(const InputParameters & par
   if (_freq_start >= _freq_end)
     mooseError("Error in " + name() +
                ". Starting frequency must be less than the ending frequency.");
-
+  // Check that frequencies are positive
+  if (_freq_start <= 0.0)
+    mooseError("Error in " + name() + ". Start and end frequencies must be positive.");
   // Check for damping
   if (_xi <= 0)
     mooseError("Error in " + name() + ". Damping ratio must be positive.");
@@ -86,6 +89,7 @@ void
 ResponseSpectraCalculator::initialize()
 {
   _frequency.clear();
+  _period.clear();
   for (VectorPostprocessorValue * ptr : _spectrum)
     ptr->clear();
 }
@@ -105,8 +109,9 @@ ResponseSpectraCalculator::execute()
     std::vector<std::vector<Real>> var_spectrum = MastodonUtils::responseSpectrum(
         _freq_start, _freq_end, _freq_num, reg_vector[1], _xi, _reg_dt);
     _frequency = var_spectrum[0];
-    *_spectrum[3 * i] = var_spectrum[1];
-    *_spectrum[3 * i + 1] = var_spectrum[2];
-    *_spectrum[3 * i + 2] = var_spectrum[3];
+    _period = var_spectrum[1];
+    *_spectrum[3 * i] = var_spectrum[2];
+    *_spectrum[3 * i + 1] = var_spectrum[3];
+    *_spectrum[3 * i + 2] = var_spectrum[4];
   }
 }