Line data Source code
1 : // MOOSE includes
2 : #include "Fragility.h"
3 : #include "PostprocessorInterface.h"
4 : #include "VectorPostprocessorInterface.h"
5 : #include "MastodonUtils.h"
6 : #include "LinearInterpolation.h"
7 : #include "DelimitedFileReader.h"
8 : #include <string>
9 :
10 : registerMooseObject("MastodonApp", Fragility);
11 :
12 : InputParameters
13 18 : Fragility::validParams()
14 : {
15 18 : InputParameters params = GeneralVectorPostprocessor::validParams();
16 36 : params.addParam<std::string>("master_file", "Name of the master file without extension.");
17 36 : params.addParam<std::string>("hazard_multiapp",
18 : "Name of the multiapp corresponding to ground motion sampling.");
19 36 : params.addParam<std::string>(
20 : "probabilistic_multiapp",
21 : "Name of the multiapp corresponding to the probabilistic simulations.");
22 36 : params.addRequiredParam<unsigned int>("num_gms",
23 : "Number of ground motions used in each intensity bin.");
24 36 : params.addParam<std::string>(
25 : "demand_variable",
26 : "Demand variable for the SSC that is also column name in the output csv "
27 : "file. Acceleration variable only."); // TODO: generalize class for all
28 : // demand variables by processing
29 : // VectorPostprocessor data and
30 : // getting rid of response spectra
31 : // calculation in this class.
32 36 : params.addParam<Real>("ssc_frequency",
33 : "Frequency at which the spectral demand of the SSC is calculated.");
34 36 : params.addParam<Real>("ssc_damping_ratio",
35 : "Damping ratio at which the spectral demand of the SSC is calculated.");
36 36 : params.addParam<Real>("dtsim", "dt for response spectra calculation.");
37 36 : params.addParam<std::string>("demand_filename",
38 : "File name that contains stochastic demand matrix. Has m x n values "
39 : "where m is number of GMs in each bin and n is num bins.");
40 36 : params.addRequiredParam<Real>(
41 : "median_capacity",
42 : "Median capacity of the SSC in terms of local demand at the SSC location.");
43 36 : params.addRequiredParam<Real>("beta_capacity", "Uncertainty in the capacity of the SSC.");
44 36 : params.addRequiredParam<unsigned int>(
45 : "num_bins", "Number of bins in the hazard curve where the risk calculation is performed.");
46 36 : params.addRequiredParam<unsigned int>("num_samples",
47 : "Number of probabilistic simulations for each bin.");
48 36 : params.addParam<unsigned int>(
49 : "num_collapses",
50 36 : 500,
51 : "Number of collapses required to calculate likelihood when using Baker's MLE.");
52 36 : params.addRequiredParam<std::vector<Real>>(
53 : "im_values", "IM values used in the bins."); // TODO: remove this later by transferring IM
54 : // values from hazard curve UO
55 36 : params.addRequiredParam<std::vector<Real>>("median_fragility_limits",
56 : "Limits for median fragility of the component.");
57 36 : params.addRequiredParam<std::vector<Real>>(
58 : "beta_fragility_limits",
59 : "Limits for the lognormal standard deviation of the component fragility.");
60 36 : params.addParam<bool>("brute_force",
61 36 : false,
62 : "Optimization "
63 : "method for fragility fitting. The following "
64 : "methods are available: brute force "
65 : "or Randomized Gradient Descent.");
66 36 : params.addParam<Real>("rgd_tolerance",
67 36 : 1e-03,
68 : "Tolerance for declaring convergence of the Randomized Gradient Descent "
69 : "algorithm.");
70 36 : params.addParam<Real>("rgd_gamma",
71 36 : 0.001,
72 : "Parameter controlling the step size of the Randomized Gradient Descent "
73 : "algorithm.");
74 36 : params.addParam<Real>("rgd_numrnd",
75 36 : 1000,
76 : "Number of random initializations in the Randomized Gradient Descent "
77 : "algorithm.");
78 36 : params.addParam<Real>("rgd_seed",
79 36 : 1028,
80 : "Seed for random number generator in the Randomized Gradient Descent "
81 : "algorithm.");
82 18 : params.addClassDescription("Calculate the seismic fragility of an SSC by postprocessing the "
83 : "results of a probabilistic or stochastic simulation.");
84 18 : return params;
85 0 : }
86 :
87 9 : Fragility::Fragility(const InputParameters & parameters)
88 : : GeneralVectorPostprocessor(parameters),
89 21 : _master_file(isParamValid("master_file") ? &getParam<std::string>("master_file") : NULL),
90 30 : _hazard_multiapp(isParamValid("hazard_multiapp") ? &getParam<std::string>("hazard_multiapp")
91 : : NULL),
92 9 : _probabilistic_multiapp(isParamValid("probabilistic_multiapp")
93 21 : ? &getParam<std::string>("probabilistic_multiapp")
94 : : NULL),
95 30 : _demand_variable(isParamValid("demand_variable") ? &getParam<std::string>("demand_variable")
96 : : NULL),
97 30 : _ssc_freq(isParamValid("ssc_frequency") ? &getParam<Real>("ssc_frequency") : NULL),
98 30 : _ssc_xi(isParamValid("ssc_damping_ratio") ? &getParam<Real>("ssc_damping_ratio") : NULL),
99 30 : _dtsim(isParamValid("dtsim") ? &getParam<Real>("dtsim") : NULL),
100 6 : _rh_file_exist(_master_file && _hazard_multiapp && _probabilistic_multiapp &&
101 15 : _demand_variable && _ssc_freq && _ssc_xi && _dtsim),
102 24 : _demand_filename(isParamValid("demand_filename") ? &getParam<std::string>("demand_filename")
103 : : NULL),
104 9 : _sd_file_exist(!!_demand_filename),
105 18 : _num_gms(getParam<unsigned int>("num_gms")),
106 18 : _median_cap(getParam<Real>("median_capacity")),
107 18 : _beta_ssc_cap(getParam<Real>("beta_capacity")),
108 18 : _num_bins(getParam<unsigned int>("num_bins")),
109 18 : _num_samples(getParam<unsigned int>("num_samples")),
110 18 : _num_collapses(getParam<unsigned int>("num_collapses")),
111 18 : _im_values(getParam<std::vector<Real>>("im_values")),
112 18 : _median_fragility_limits(getParam<std::vector<Real>>("median_fragility_limits")),
113 18 : _beta_fragility_limits(getParam<std::vector<Real>>("beta_fragility_limits")),
114 9 : _im(declareVector("intensity")),
115 9 : _median_demand(declareVector("demand_median")),
116 9 : _beta_demand(declareVector("demand_beta")),
117 9 : _conditional_pf(declareVector("conditional_pf")),
118 9 : _median_fragility(declareVector("fragility_median")),
119 9 : _beta_fragility(declareVector("fragility_beta")),
120 9 : _loglikelihood(declareVector("loglikelihood")),
121 18 : _brute_force(getParam<bool>("brute_force")),
122 18 : _rgd_tolerance(getParam<Real>("rgd_tolerance")),
123 18 : _rgd_gamma(getParam<Real>("rgd_gamma")),
124 18 : _rgd_numrnd(getParam<Real>("rgd_numrnd")),
125 27 : _rgd_seed(getParam<Real>("rgd_seed"))
126 : {
127 9 : if (_rh_file_exist && _sd_file_exist)
128 0 : mooseError("Error in block '" + name() +
129 : "'. Both response history file options and stochastic demand file options "
130 : "are provided. Provide exactly one set of them.");
131 9 : if (!_rh_file_exist && !_sd_file_exist)
132 0 : mooseError(
133 0 : "Error in block '" + name() +
134 : "'. Either one or more of the response history file options are missing or "
135 : "none of the response history file options or stochastic demand file "
136 : "options are provided. Provide exactly one of them. \n\n **If response history files "
137 : "are to be used, please provide the input parameters, master_file, hazard_multiapp, "
138 : "probabilistic_multiapp, demand_variable, ssc_freq, ssc_xi, and dt.");
139 :
140 9 : if (!_rh_file_exist && (_master_file || _hazard_multiapp || _probabilistic_multiapp ||
141 3 : _demand_variable || _ssc_freq || _ssc_xi || _dtsim))
142 0 : mooseError(
143 0 : "Error in block '" + name() +
144 : "'. If response history files "
145 : "are to be used, please provide ALL of the input parameters, master_file, hazard_multiapp, "
146 : "probabilistic_multiapp, demand_variable, ssc_freq, ssc_xi, and dt.");
147 :
148 9 : if (_rh_file_exist)
149 : {
150 : // Check for non-positive SSC frequency
151 6 : if (*_ssc_freq <= 0)
152 0 : mooseError("Error in block '" + name() + "'. SSC frequency must be positive.");
153 : // Check for non-positive SSC damping
154 6 : if (*_ssc_xi <= 0)
155 0 : mooseError("Error in block '" + name() + "'. SSC damping ratio must be positive.");
156 : }
157 :
158 : // // Check for non-positive median and beta of capacity
159 9 : if (_median_cap <= 0 || _beta_ssc_cap <= 0)
160 0 : mooseError(
161 0 : "Error in block '" + name() +
162 : "'. Median and lognormal standard deviation of the capacity of the SSC must be positive.");
163 : // Check if median and lognormal standard deviation limits for fraglity are of size 2
164 9 : if (_median_fragility_limits.size() != 2 || _beta_fragility_limits.size() != 2)
165 0 : mooseError("Error in block '" + name() +
166 : "'. Please provide two limits (minimum and maximum) "
167 : "for median and lognormal standard deviation of the "
168 : "component.");
169 9 : if (_median_fragility_limits[0] >= _median_fragility_limits[1])
170 0 : mooseError("Error in block '" + name() +
171 : "'. Minimum median fragility should be less than the maximum median fragility.");
172 9 : if (_beta_fragility_limits[0] >= _beta_fragility_limits[1])
173 0 : mooseError(
174 0 : "Error in block '" + name() +
175 : "'. Minimum should be less than the maximum for fragility lognormal standard deviation.");
176 9 : if (_median_fragility_limits[0] <= 0 || _beta_fragility_limits[0] <= 0)
177 0 : mooseError("Error in block '" + name() +
178 : "'. Limits of median and beta of the component fragility must be positive.");
179 : // Check if _im_values is of the same size as the number of bins
180 9 : if (_im_values.size() != _num_bins)
181 0 : mooseError("Error in block '" + name() +
182 : "'. Number of IM values should be the same as the number of bins.");
183 9 : }
184 :
185 : void
186 9 : Fragility::initialize()
187 : {
188 9 : _median_demand.clear();
189 9 : _beta_demand.clear();
190 9 : _im.clear();
191 9 : _conditional_pf.clear();
192 9 : _median_fragility.clear();
193 9 : _beta_fragility.clear();
194 9 : _loglikelihood.clear();
195 9 : }
196 :
197 : void
198 9 : Fragility::execute()
199 : {
200 9 : _im.resize(_num_bins);
201 9 : _median_demand.resize(_num_bins);
202 9 : _beta_demand.resize(_num_bins);
203 9 : _conditional_pf.resize(_num_bins);
204 9 : _median_fragility.resize(1);
205 9 : _beta_fragility.resize(1);
206 9 : _loglikelihood.resize(1);
207 51 : for (std::size_t bin = 0; bin < _num_bins; bin++)
208 : {
209 42 : _im[bin] = _im_values[bin];
210 : std::vector<Real> stoc_demands;
211 42 : if (_sd_file_exist)
212 : // reading demands from stochastic demand csv files if it is a restart PRA
213 36 : stoc_demands = Fragility::readDemandsFromSDFiles(bin);
214 : else
215 : // calculating demands from RH csv files if it is a restart PRA
216 48 : stoc_demands = Fragility::calcDemandsFromRHFiles(bin);
217 42 : _median_demand[bin] = MastodonUtils::median(stoc_demands);
218 42 : _beta_demand[bin] = MastodonUtils::lognormalStandardDeviation(stoc_demands);
219 :
220 42 : Real demand_median = _median_demand[bin];
221 : Real demand_scale = _beta_demand[bin];
222 42 : Real capacity_median = _median_cap;
223 42 : Real capacity_scale = _beta_ssc_cap;
224 :
225 42 : _conditional_pf[bin] = MastodonUtils::greaterProbability(
226 : demand_median, demand_scale, capacity_median, capacity_scale);
227 42 : _console << "**********\nFinished calculations for bin: " << bin
228 42 : << " \n Median demand: " << _median_demand[bin]
229 42 : << " \n Lognormal standard deviation: " << _beta_demand[bin]
230 42 : << " \n Conditional probability of failure: " << _conditional_pf[bin]
231 42 : << "\n**********\n";
232 42 : }
233 9 : std::vector<Real> fitted_vals = MastodonUtils::maximizeLogLikelihood(_im,
234 9 : _conditional_pf,
235 : _median_fragility_limits,
236 : _beta_fragility_limits,
237 : _num_collapses,
238 : _brute_force,
239 9 : _rgd_tolerance,
240 9 : _rgd_gamma,
241 9 : _rgd_numrnd,
242 9 : _rgd_seed);
243 9 : _median_fragility[0] = fitted_vals[0];
244 9 : _beta_fragility[0] = fitted_vals[1];
245 9 : _loglikelihood[0] = MastodonUtils::calcLogLikelihood(
246 9 : _im, _conditional_pf, _median_fragility[0], _beta_fragility[0], _num_collapses);
247 9 : }
248 :
249 : std::vector<Real>
250 18 : Fragility::readDemandsFromSDFiles(unsigned int bin)
251 : {
252 18 : _console << "Stochastic demands file exists. Reading from " << *_demand_filename << "\n";
253 : // Reading stochastic demand vector from the stochastic demands file of
254 : // the component. The file has m x n demand values, where m is the
255 : // number of GMs per bin and n is the number of bins.
256 : std::vector<Real> stoc_demands;
257 18 : stoc_demands.resize(_num_samples * _num_gms);
258 18 : MooseUtils::DelimitedFileReader stoc_demands_file(*_demand_filename);
259 18 : stoc_demands_file.read();
260 18 : if ((stoc_demands_file.getNames()).size() != _num_bins)
261 0 : mooseError("Error in block '" + name() +
262 : "'. Number of columns in stochastic demands file is not the "
263 : "same as the number of bins.");
264 36 : stoc_demands = stoc_demands_file.getData("bin_" + std::to_string(bin + 1));
265 18 : if (stoc_demands.size() != _num_gms * _num_samples)
266 0 : mooseError("Error in block '" + name() +
267 : "'. Number of rows in stochastic demands file is not the "
268 : "same as the product of the number of GMs and number of samples.");
269 18 : return stoc_demands;
270 18 : }
271 :
272 : std::vector<Real>
273 24 : Fragility::calcDemandsFromRHFiles(unsigned int bin)
274 : {
275 : std::vector<Real> demand_sample;
276 : std::vector<std::vector<Real>> demand_sample_spectrum;
277 : std::vector<Real> stoc_demands;
278 24 : stoc_demands.resize(_num_samples * _num_gms);
279 : std::vector<Real> demand_time;
280 : std::string demand_sample_filename;
281 : Real k = 0;
282 96 : for (unsigned int i = 0; i < _num_gms; ++i)
283 : {
284 288 : for (unsigned int j = 0; j < _num_samples; ++j)
285 : {
286 432 : demand_sample_filename = *_master_file + "_out_" + *_hazard_multiapp +
287 432 : MastodonUtils::zeropad(bin * _num_gms + i, _num_bins * _num_gms) +
288 864 : "_" + *_probabilistic_multiapp + std::to_string(j) + ".csv";
289 432 : _console << "In block '" + name() + "'. Reading file: " << demand_sample_filename
290 216 : << std::endl;
291 216 : MooseUtils::DelimitedFileReader demand_sample_file(demand_sample_filename);
292 216 : demand_sample_file.read();
293 216 : demand_sample = demand_sample_file.getData(*_demand_variable);
294 216 : demand_time = demand_sample_file.getData("time");
295 216 : demand_sample = MastodonUtils::regularize(
296 216 : demand_sample, demand_time, *_dtsim)[1]; // regularize the demand sample
297 : demand_sample_spectrum =
298 432 : MastodonUtils::responseSpectrum(0.01, 100, 401, demand_sample, *_ssc_xi, *_dtsim);
299 216 : LinearInterpolation spectraldemand(demand_sample_spectrum[0], demand_sample_spectrum[4]);
300 216 : stoc_demands[k] = spectraldemand.sample(*_ssc_freq);
301 216 : k++;
302 216 : }
303 : }
304 24 : return stoc_demands;
305 24 : }
306 :
307 : // TODO: Currently all stochastic simulations have the same termination time.
308 : // Add capability to terminate at various times.
309 : // TODO: Add capability to calculate demands from transferring VectorPostprocessor
310 : // data instead of reading csv outputs.
|
