Failure analysis

Computational Model Description

This section outlines the setup and execution of a Parallel Subset Simulation (PSS) framework using a 2D MSRE (Molten Salt Reactor Experiment) geometry. The framework utilizes two input files: the main PSS input file and the finite element simulation input file. The PSS input file generates combinations of input parameters and passes them to the simulation file, which performs the finite element analysis to determine the maximum stress corresponding to a user-defined infiltration amount. The PSS input file then retrieves the simulation results, such as the maximum stress, for further analysis. It is important to note that the finite element simulation input file adheres to a setup analogous to the one described in the 3D stress analysis, but is specifically configured for 2D simulations. Consequently, the detailed discussion of the finite element setup is omitted here for brevity.

Files used by this model include:

• MOOSE input files: pss.i and Baseline_Input.i

• Exodus reference solution file for initializing the infiltration amount

• Mesh input file

• A python code for postprocessing

• HPC job file to submit the simulation

This document reviews the important elements of the input file (pss.i), listed in full here:

# ==============================================================================
# Parellel subset simulation setup to perform rare-event analysis
# Application : MOOSE
# ------------------------------------------------------------------------------
# Idaho Falls, INL, 2025
# Author(s): V Prithivirajan, Ben Spencer
# If using or referring to this model, please cite as explained on
# https://mooseframework.inl.gov/virtual_test_bed/citing.html
# ==============================================================================

[StochasticTools] #Configures the overall stochastic analysis and settings
[]

[Distributions] # Defines distributions for input paramters
  [INF] #Infiltration amount
    type = Uniform
    lower_bound = 0
    upper_bound = 1
  []

  [E] #Youngs modulus [Pa]
    type = Uniform
    lower_bound = 9e9
    upper_bound = 15e9
  []
  [K] #Thermal conductivity [W/mK]
    type = Uniform
    lower_bound = 25
    upper_bound = 100
  []
  [PD] #Power density [W/m^3]
    type = Uniform
    lower_bound = 2e6
    upper_bound = 5e7
  []
  [nu] #Poisson ratio
    type = Uniform
    lower_bound = 0.13
    upper_bound = 0.21
  []
  [Tinf] #Salt temperature [K]
    type = Uniform
    lower_bound = 823
    upper_bound = 1023
  []
  [htc] #Heat transfer coefficient [W/m^2K]
    type = Uniform
    lower_bound = 3500
    upper_bound = 5500
  []
  [CTE] #Coefficient of thermal expansion [1/K]
    type = Uniform
    lower_bound = 3.5e-6
    upper_bound = 6.0e-6
  []
[]

[Samplers] # Defines the sampler methods
  [sample]
    type = ParallelSubsetSimulation
    distributions = 'INF E K PD nu Tinf htc CTE'
    num_samplessub = 1000 #Number of simulations within a subset
    num_subsets = 5 #total number of subsets
    num_parallel_chains = 100 #Number of parellel simulations run concurrently
    output_reporter = 'constant/reporter_transfer:maxstress:value'
    inputs_reporter = 'adaptive_MC/inputs'
    seed = 1012
    execute_on = 'PRE_MULTIAPP_SETUP'
  []
[]

[MultiApps] #Manages the execution of sub-applications within the main simulation
  [sub]
    type = SamplerFullSolveMultiApp
    input_files = 'Baseline_Input.i'
    sampler = sample
  []
[]

[Transfers] # Facilitates transfer of data between main and sub app
  [reporter_transfer]
    type = SamplerReporterTransfer
    from_reporter = 'maxstress/value'
    stochastic_reporter = 'constant'
    from_multi_app = sub
    sampler = sample
  []
[]

[Controls]
  [cmdline] #Passes inputs to simulation based on sampler output
    type = MultiAppSamplerControl
    multi_app = sub
    sampler = sample
    param_names = 'INF E K PD nu Tinf htc CTE'
  []
[]

[Reporters]
  [constant]
    type = StochasticReporter
    outputs = none
  []
  [adaptive_MC] # Guides input selection towards rare-event
    type = AdaptiveMonteCarloDecision
    output_value = constant/reporter_transfer:maxstress:value
    inputs = 'inputs'
    sampler = sample
  []
[]

[Executioner]
  type = Transient
[]

[Outputs]
  [out]
    type = JSON
  []
[]

[StochasticTools]
  #Configures the overall stochastic analysis and settings
[]

[Distributions]
  # Defines distributions for input paramters
  [INF]
    #Infiltration amount
    type = Uniform
    lower_bound = 0
    upper_bound = 1
  []

  [E]
    #Youngs modulus [Pa]
    type = Uniform
    lower_bound = 9e9
    upper_bound = 15e9
  []
  [K]
    #Thermal conductivity [W/mK]
    type = Uniform
    lower_bound = 25
    upper_bound = 100
  []
  [PD]
    #Power density [W/m^3]
    type = Uniform
    lower_bound = 2e6
    upper_bound = 5e7
  []
  [nu]
    #Poisson ratio
    type = Uniform
    lower_bound = 0.13
    upper_bound = 0.21
  []
  [Tinf]
    #Salt temperature [K]
    type = Uniform
    lower_bound = 823
    upper_bound = 1023
  []
  [htc]
    #Heat transfer coefficient [W/m^2K]
    type = Uniform
    lower_bound = 3500
    upper_bound = 5500
  []
  [CTE]
    #Coefficient of thermal expansion [1/K]
    type = Uniform
    lower_bound = 3.5e-6
    upper_bound = 6.0e-6
  []
[]

[Samplers]
  # Defines the sampler methods
  [sample]
    type = ParallelSubsetSimulation
    distributions = 'INF E K PD nu Tinf htc CTE'
    num_samplessub = 1000 #Number of simulations within a subset
    num_subsets = 5 #total number of subsets
    num_parallel_chains = 100 #Number of parellel simulations run concurrently
    output_reporter = 'constant/reporter_transfer:maxstress:value'
    inputs_reporter = 'adaptive_MC/inputs'
    seed = 1012
    execute_on = 'PRE_MULTIAPP_SETUP'
  []
[]

[MultiApps]
  #Manages the execution of sub-applications within the main simulation
  [sub]
    type = SamplerFullSolveMultiApp
    input_files = 'Baseline_Input.i'
    sampler = sample
  []
[]

[Transfers]
  # Facilitates transfer of data between main and sub app
  [reporter_transfer]
    type = SamplerReporterTransfer
    from_reporter = 'maxstress/value'
    stochastic_reporter = 'constant'
    from_multi_app = sub
    sampler = sample
  []
[]

[Controls]
  [cmdline]
    #Passes inputs to simulation based on sampler output
    type = MultiAppSamplerControl
    multi_app = sub
    sampler = sample
    param_names = 'INF E K PD nu Tinf htc CTE'
  []
[]

[Reporters]
  [constant]
    type = StochasticReporter
    outputs = none
  []
  [adaptive_MC]
    # Guides input selection towards rare-event
    type = AdaptiveMonteCarloDecision
    output_value = constant/reporter_transfer:maxstress:value
    inputs = 'inputs'
    sampler = sample
  []
[]

Running the model

This analysis requires HPC resources. To run this model using the MOOSE combined module executable, run the job file in a HPC system:

#!/bin/bash

#PBS -N PSS
#PBS -l select=25:ncpus=48:mpiprocs=40
#PBS -l walltime=05:00:00
#PBS -P nrc

# Change to the directory from which the job was submitted
cd $PBS_O_WORKDIR

# Load the necessary modules
module purge
module load use.moose moose-dev

mpiexec -n 1000 path_to_combined-opt -i pss.i | tee progress_log.txt 

Note: Please note that this script was written for a computing cluster using the PBS scheduling system. Therefore, the job script is based on PBS commands. If you need to use SLURM, ensure to adjust the script accordingly.

The following output files will be produced:

• The JSON results file: pss_out.json

This is analyzed using the following python file:

import json
import numpy as np
import os
import seaborn as sns
import matplotlib.pyplot as plt
import random
from scipy.stats import entropy

cwd = os.getcwd()

# Define number of bins for the two histograms
bins1 = np.array([20, 20, 20, 20, 20, 20, 20, 20])  # Number of bins for Monte Carlo samples
bins2 = np.array([4, 10, 4, 5, 15, 15, 15, 6])  # Number of bins for failed samples

os.chdir(cwd)

# Define file path
file_path = cwd + "/pss_out.json"

Nsubsets_to_consider = 7000

# Initialize lists to store inputs and outputs
all_inputs = []
all_outputs = []

# Load the JSON file
with open(file_path, "r") as f:
    data = json.load(f)

# Extract data for each time step
time_steps = data["time_steps"]  # Adjust key if necessary
length_list = len(time_steps)

for i in range(1, length_list):
    step = time_steps[i]
    inputs = step["adaptive_MC"]["inputs"]  # Access the inputs
    transposed_inputs = np.array(inputs).T
    output = step["adaptive_MC"]["output_required"]  # Access the output

    # Append inputs and outputs
    all_inputs.extend(transposed_inputs)
    all_outputs.extend(output)

# Convert to NumPy arrays
inputs_array = np.array(all_inputs).reshape(-1, 8)
outputs_array = np.array(all_outputs)
inputs_array = inputs_array[:Nsubsets_to_consider,:]
outputs_array = outputs_array[:Nsubsets_to_consider]

# Check shapes
print("Inputs shape:", inputs_array.shape)
print("Outputs shape:", outputs_array.shape)

# Failure probability
last_1000_points = outputs_array[-1000:]
subset_prob_calc = len(last_1000_points[last_1000_points > 15e6]) / len(last_1000_points)
global_prob = subset_prob_calc * 1e-6
ind1 = outputs_array > 15e6
failed_points = len(outputs_array[ind1])

plt.figure(dpi=300)
plt.plot(outputs_array * 1e-6)
plt.xlabel('No. of samples')
plt.ylabel('Max. Stress (MPa)')
plt.grid(False)
plt.show()

rndvals = random.sample(range(1000), failed_points)
kl_divergence = np.zeros(8)

# Use LaTeX for all text rendering
# plt.rcParams['text.usetex'] = True

labels = [
    'Infiltration (\%)',
    r"Young's Modulus (GPa)",
    r'Thermal Conductivity (W/mK)',
    r'Power Density (MW/m$^3$)',
    r'Poisson Ratio',
    r'Interface Temperature (K)',
    r'Heat Transfer Coeff. (W/m$^2$K)',
    r'Coeff. of Thermal Expansion (1/K)'
]

bounds = [
    (0, 1),
    (9, 15),
    (25, 100),
    (10, 250),
    (0.13, 0.21),
    (823, 1023),
    (3500, 5500),
    (3.5e-6, 6.0e-6)
]

for idx0 in range(8):
    hist1 = inputs_array[rndvals, idx0]
    hist2 = inputs_array[ind1, idx0]

    # Normalize the histograms to create probability distributions
    prob_dist1 = hist1 / np.sum(hist1)
    prob_dist2 = hist2 / np.sum(hist2)

    if idx0 == 3:
        fact = 1 / (0.2 * 1e6)
    elif idx0 == 1:
        fact = 1e-9
    else:
        fact = 1

    sns.set(style="whitegrid")
    plt.figure(figsize=(4, 4), dpi=300)

    plt.hist(inputs_array[0:1000, idx0] * fact, bins=bins1[idx0], density=True, alpha=0.5, color="b", label="Monte Carlo Samples")
    plt.hist(inputs_array[ind1, idx0] * fact, bins=bins2[idx0], density=True, alpha=0.5, color="r", label="Failed Samples")

    plt.xlabel(labels[idx0])
    plt.ylabel("Density")
    plt.legend()
    plt.grid(False)
    plt.show()

# Extract the input and output data for the failed points
failed_inputs = inputs_array[ind1]
failed_outputs = outputs_array[ind1]

# Define the multiplication factors
factors = np.ones(8)
factors[1] = 1e-9
factors[3] = 1 / (0.2 * 1e6)

# Apply the multiplication factors to the inputs
adjusted_inputs = failed_inputs * factors

# Apply the multiplication factor to the output
adjusted_outputs = failed_outputs * 1e-6

# Combine the adjusted input and output data into one matrix
combined_matrix = np.hstack((adjusted_inputs, adjusted_outputs.reshape(-1, 1)))

# Initialize arrays to store the minimum and maximum values
min_values = np.zeros(8)
max_values = np.zeros(8)

# Calculate the minimum and maximum values for each column
for idx in range(8):
    min_values[idx] = np.min(combined_matrix[:, idx])
    max_values[idx] = np.max(combined_matrix[:, idx])

# Display the minimum and maximum values
print("Min values:")
for idx, min_val in enumerate(min_values):
    if idx == 7:
        print(f"{min_val:.4e}")
    else:
        print(f"{min_val:.4f}")

print("\nMax values:")
for idx, max_val in enumerate(max_values):
    if idx == 7:
        print(f"{max_val:.4e}")
    else:
        print(f"{max_val:.4f}")

This Python script reads a JSON output file and generates histograms of input variables corresponding to failure, defined by a maximum stress exceeding 15 MPa. It compares these histograms to those from random Monte Carlo samples for each input variable. Additionally, the script provides the minimum and maximum ranges of input values associated with this failure.