# Using Stochastic Tools with Multiphysics Models

The purpose of this document is to present a multiphysics example using the Stochastic Tools Module. The intention is to showcase the capabilities of the module to produce statistically relevant results including uncertainty propagation and sensitivity, as well as the module's surrogate modeling infrastructure. The problem of interest is a thermomechanics model using a combination of the Heat Conduction and Tensor Mechanics modules. The problem involves multiple uncertain material properties and multiple quantities of interest (QoI). Using both Monte Carlo sampling and polynomial chaos surrogate modeling, the effect of these properties' uncertainties are quantified with uncertainty propagation and global sensitivity analysis.

## Problem Description

The problem of interest involves a steady-state thermomechanics model. The geometry is a 3-D finite hollow cylinder with two concentric layers of different material properties, seen in Figure 1. Due to symmetry, only 1/8 of the cylinder is represented by the mesh. The inner surface of the ring is exposed to a surface heat source that has a cosine shape along the axis of the cylinder with its peak being at the center. The end and outside of the cylinder have convective boundary conditions. The cylinder is free to displace in all directions due to the thermal expansion. The relevant material and geometric properties are listed in Table 1. This table also lists the "uncertain" parameters that will be described in the following section. For reference, the temperature and displacement profiles are shown in Figure 4 to Figure 8 where the uncertain properties are set to some arbitrary values. The full input file for this model is shown in Listing 1.

Figure 1: Model Problem Geometry

Table 1: Material Properties for Thermomechanics Cylinder

PropertySymbolValueUnits
Half Cyliner Height3m
Inner Ring Width0.1m
Outer Ring Width0.1m
Outer Heat Transfer Coef.10W/mK
Outer Free Temperature300K
End Heat Transfer Coef.10W/mK
End Free Temperature300K
Heat Source MagnitudeUncertainW
Inner Thermal ConductivityUncertainW/mK
Outer Thermal ConductivityUncertainW/mK
Inner Young's ModulusUncertainPa
Outer Young's ModulusUncertainPa
Inner Poisson's RatioUncertain
Outer Poisson's RatioUncertain
Inner Thermal Expansion Coef.Uncertain1/K
Outer Thermal Expansion Coef.Uncertain1/K
At Rest Temperature300K

Figure 4: Temperature (K)

Figure 5: Displacement (m)

Figure 6: X-Displacement (m)

Figure 7: Y-Displacement (m)

Figure 8: Z-Displacement (m)

Listing 1: Thermomechanics model input file

# Generate 1/4 of a 2-ring disk and extrude it by half to obtain
# 1/8 of a 3D tube.  Mirror boundary conditions will exist on the
# cut portions.

[Mesh]
[disk]
type = ConcentricCircleMeshGenerator
num_sectors = 10
rings = '1 1 1'
has_outer_square = false
preserve_volumes = false
portion = top_right
[]
[ring]
type = BlockDeletionGenerator
input = disk
block_id = 1
new_boundary = 'inner'
[]
[cylinder]
type = MeshExtruderGenerator
input = ring
extrusion_vector = '0 0 1.5'
num_layers = 15
bottom_sideset = 'back'
top_sideset = 'front'
[]
[]

[Variables]
[T]
initial_condition = 300
[]
[disp_x]
[]
[disp_y]
[]
[disp_z]
[]
[]

[Kernels]
[hc]
type = HeatConduction
variable = T
[]
[TensorMechanics]
displacements = 'disp_x disp_y disp_z'
[]
[]

[BCs]
[temp_inner]
type = FunctionNeumannBC
variable = T
boundary = 'inner'
function = surface_source
[]
[temp_front]
type = ConvectiveHeatFluxBC
variable = T
boundary = 'front'
T_infinity = 300
heat_transfer_coefficient = 10
[]
[temp_outer]
type = ConvectiveHeatFluxBC
variable = T
boundary = 'outer'
T_infinity = 300
heat_transfer_coefficient = 10
[]
# mirror boundary conditions.
[disp_x]
type = DirichletBC
variable = disp_x
boundary = 'left'
value = 0.0
[]
[disp_y]
type = DirichletBC
variable = disp_y
boundary = 'bottom'
value = 0.0
[]
[disp_z]
type = DirichletBC
variable = disp_z
boundary = 'back'
value = 0.0
[]
[]

[Materials]
[cond_inner]
type = GenericConstantMaterial
block = 2
prop_names = thermal_conductivity
prop_values = 25
[]
[cond_outer]
type = GenericConstantMaterial
block = 3
prop_names = thermal_conductivity
prop_values = 100
[]
[elasticity_tensor_inner]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 2.1e5
poissons_ratio = 0.3
block = 2
[]
[elasticity_tensor_outer]
type = ComputeIsotropicElasticityTensor
youngs_modulus = 3.1e5
poissons_ratio = 0.2
block = 3
[]
[thermal_strain_inner]
type = ComputeThermalExpansionEigenstrain
thermal_expansion_coeff = 2e-6
temperature = T
stress_free_temperature = 300
eigenstrain_name = eigenstrain_inner
block = 2
[]
[thermal_strain_outer]
type = ComputeThermalExpansionEigenstrain
thermal_expansion_coeff = 1e-6
temperature = T
stress_free_temperature = 300
eigenstrain_name = eigenstrain_outer
block = 3
[]
[strain_inner] #We use small deformation mechanics
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = 'eigenstrain_inner'
block = 2
[]
[strain_outer] #We use small deformation mechanics
type = ComputeSmallStrain
displacements = 'disp_x disp_y disp_z'
eigenstrain_names = 'eigenstrain_outer'
block = 3
[]
[stress] #We use linear elasticity
type = ComputeLinearElasticStress
[]
[]

[Functions]
[surface_source]
type = ParsedFunction
value = 'Q_t*pi/2.0/3.0 * cos(pi/3.0*z)'
vars = 'Q_t'
vals = heat_source
[]
[]

[Executioner]

#Preconditioned JFNK (default)
solve_type = 'PJFNK'
petsc_options_iname = '-pc_type -pc_hypre_type -ksp_gmres_restart'
petsc_options_value = 'hypre boomeramg 101'

l_max_its = 30
nl_max_its = 100
nl_abs_tol = 1e-9
l_tol = 1e-04
[]

[Preconditioning]
[SMP]
type = SMP
full = true
[]
[]

[Controls]
[stochastic]
[]
[]

[VectorPostprocessors]
[temp_center]
type = LineValueSampler
variable = T
start_point = '1 0 0'
end_point = '1.2 0 0'
num_points = 11
sort_by = 'x'
[]
[temp_end]
type = LineValueSampler
variable = T
start_point = '1 0 1.5'
end_point = '1.2 0 1.5'
num_points = 11
sort_by = 'x'
[]
[dispx_center]
type = LineValueSampler
variable = disp_x
start_point = '1 0 0'
end_point = '1.2 0 0'
num_points = 11
sort_by = 'x'
[]
[dispx_end]
type = LineValueSampler
variable = disp_x
start_point = '1 0 1.5'
end_point = '1.2 0 1.5'
num_points = 11
sort_by = 'x'
[]
[dispz_end]
type = LineValueSampler
variable = disp_z
start_point = '1 0 1.5'
end_point = '1.2 0 1.5'
num_points = 11
sort_by = 'x'
[]
[]

[Postprocessors]
[heat_source]
type = FunctionValuePostprocessor
function = 1
scale_factor = 10000
execute_on = linear
[]
[temp_center_inner]
type = PointValue
variable = T
point = '1 0 0'
[]
[temp_center_outer]
type = PointValue
variable = T
point = '1.2 0 0'
[]
[temp_end_inner]
type = PointValue
variable = T
point = '1 0 1.5'
[]
[temp_end_outer]
type = PointValue
variable = T
point = '1.2 0 1.5'
[]
[dispx_center_inner]
type = PointValue
variable = disp_x
point = '1 0 0'
[]
[dispx_center_outer]
type = PointValue
variable = disp_x
point = '1.2 0 0'
[]
[dispx_end_inner]
type = PointValue
variable = disp_x
point = '1 0 1.5'
[]
[dispx_end_outer]
type = PointValue
variable = disp_x
point = '1.2 0 1.5'
[]
[dispz_inner]
type = PointValue
variable = disp_z
point = '1 0 1.5'
[]
[dispz_outer]
type = PointValue
variable = disp_z
point = '1.2 0 1.5'
[]
[]

[Outputs]
exodus = false
csv = false
[]

(modules/combined/examples/stochastic/graphite_ring_thermomechanics.i)

### Uncertain Parameters

A total of nine properties of the model are not known exactly, but with some known probability of values occurring. The probabilities are represented by each parameter's probability density function. All parameters have a independent uniform distribution, , where and are the lower and upper limit values for the property, respectively. Table 2 lists the details of each of the property's distribution.

Table 2: Uniform distribution parameters for uncertain properties

IndexProperty
12030
290110
3900011000
41.03.0
50.51.5
62.02.2
73.03.2
80.290.31
90.190.21

### Quantities of Interest

There are a total of ten QoIs for the model, which involve temperature and displacement:

1. Temperature at center of inner surface —

2. Temperature at center of outer surface —

3. Temperature at end of inner surface —

4. Temperature at end of outer surface —

5. x-displacement at center of inner surface —

6. x-displacement at center of outer surface —

7. x-displacement at end of inner surface —

8. x-displacement at end of outer surface —

9. z-displacement at end of inner surface —

10. z-displacement at end of outer surface —

Figure 2 shows geometrically where these QoIs are located in the model.

Figure 2: Geometric description of quantities of interest

## Results

In this exercise, we will use the statistics and Sobol sensitivity capabilities available in the stochastic tools module. The goal of this exercise is to understand how the uncertainty in the parameters affects the the resulting QoIs. This is done through sampling the model at different perturbations of the parameters and performing statistical calculations on resulting QoI values. Two methods are used to perform this analysis. First is using the sampler system to perturb the uncertain properties and retrieve the QoIs which will undergo the analysis. The second is training a polynomial chaos surrogate and using that reduced order model to sample and perform the analysis. The idea is that many evaluations of the model are necessary to compute accurate statistical quantities and surrogate modeling speeds up this computation by requiring much fewer full model evaluations for training and is significantly faster to evaluate once trained.

Using latin hypercube sampling, the thermomechanics model was run with a total of 100,000 samples, the input file is shown by Listing 2. A order four polynomial chaos surrogate was training using a Smolyak sparse quadrature for a total of 7,344 runs of the full model. The training input is shown by Listing 3 and the evaluation input is shown by Listing 4. Table 3 shows the run-time for sampling the full order model and training and evaluating the surrogate. We see here that cumulative time for training and evaluating the surrogate is much smaller than just sampling the full order model, this is because building the surrogate required far fewer evaluations of the full model and evaluating the surrogate is much faster than evaluating the full model.

Table 3: Stochastic run-time results

SimulationSamplesCPU Time
Full-Order Sampling100,000176 hr
Polynomial Chaos — Training7,34413.7 hr
Polynomial Chaos — Evaluation100,0006.8 s

Listing 2: Latin hypercube sampling and statistics input file

[StochasticTools]
[]

[Distributions]
[cond_inner]
type = Uniform
lower_bound = 20
upper_bound = 30
[]
[cond_outer]
type = Uniform
lower_bound = 90
upper_bound = 110
[]
[heat_source]
type = Uniform
lower_bound = 9000
upper_bound = 11000
[]
[alpha_inner]
type = Uniform
lower_bound = 1e-6
upper_bound = 3e-6
[]
[alpha_outer]
type = Uniform
lower_bound = 5e-7
upper_bound = 1.5e-6
[]
[ymod_inner]
type = Uniform
lower_bound = 2e5
upper_bound = 2.2e5
[]
[ymod_outer]
type = Uniform
lower_bound = 3e5
upper_bound = 3.2e5
[]
[prat_inner]
type = Uniform
lower_bound = 0.29
upper_bound = 0.31
[]
[prat_outer]
type = Uniform
lower_bound = 0.19
upper_bound = 0.21
[]
[]

[Samplers]
[sample]
type = LatinHypercube
num_rows = 100000
distributions = 'cond_inner cond_outer heat_source alpha_inner alpha_outer ymod_inner ymod_outer prat_inner prat_outer'
execute_on = INITIAL
num_bins = 10
[]
[]

[MultiApps]
[sub]
type = SamplerFullSolveMultiApp
input_files = graphite_ring_thermomechanics.i
sampler = sample
mode = batch-reset
[]
[]

[Transfers]
[sub]
type = SamplerParameterTransfer
multi_app = sub
sampler = sample
parameters = 'Materials/cond_inner/prop_values Materials/cond_outer/prop_values
Postprocessors/heat_source/scale_factor
Materials/thermal_strain_inner/thermal_expansion_coeff Materials/thermal_strain_outer/thermal_expansion_coeff
Materials/elasticity_tensor_inner/youngs_modulus Materials/elasticity_tensor_outer/youngs_modulus
Materials/elasticity_tensor_inner/poissons_ratio Materials/elasticity_tensor_outer/poissons_ratio'
to_control = 'stochastic'
check_multiapp_execute_on = false
[]
[data]
type = SamplerPostprocessorTransfer
multi_app = sub
sampler = sample
to_vector_postprocessor = storage
from_postprocessor = 'temp_center_inner  temp_center_outer  temp_end_inner  temp_end_outer
dispx_center_inner dispx_center_outer dispx_end_inner dispx_end_outer
dispz_inner dispz_outer'
[]
[]

[VectorPostprocessors]
[storage]
type = StochasticResults
[]
[stats]
type = Statistics
vectorpostprocessors = 'storage'
compute = 'mean stddev'
ci_method = 'percentile'
ci_levels = '0.05'
[]
[]

[Outputs]
csv = true
execute_on = TIMESTEP_END
[]

(modules/combined/examples/stochastic/lhs_uniform.i)

Listing 3: Polynomial chaos training input file

[StochasticTools]
[]

[Distributions]
[cond_inner]
type = Uniform
lower_bound = 20
upper_bound = 30
[]
[cond_outer]
type = Uniform
lower_bound = 90
upper_bound = 110
[]
[heat_source]
type = Uniform
lower_bound = 9000
upper_bound = 11000
[]
[alpha_inner]
type = Uniform
lower_bound = 1e-6
upper_bound = 3e-6
[]
[alpha_outer]
type = Uniform
lower_bound = 5e-7
upper_bound = 1.5e-6
[]
[ymod_inner]
type = Uniform
lower_bound = 2e5
upper_bound = 2.2e5
[]
[ymod_outer]
type = Uniform
lower_bound = 3e5
upper_bound = 3.2e5
[]
[prat_inner]
type = Uniform
lower_bound = 0.29
upper_bound = 0.31
[]
[prat_outer]
type = Uniform
lower_bound = 0.19
upper_bound = 0.21
[]
[]

[GlobalParams]
distributions = 'cond_inner cond_outer heat_source alpha_inner alpha_outer ymod_inner ymod_outer prat_inner prat_outer'
[]

[Samplers]
[sample]
sparse_grid = smolyak
order = 5
execute_on = INITIAL
[]
[]

[MultiApps]
[sub]
type = SamplerFullSolveMultiApp
input_files = graphite_ring_thermomechanics.i
sampler = sample
mode = batch-reset
[]
[]

[Transfers]
[sub]
type = SamplerParameterTransfer
multi_app = sub
sampler = sample
parameters = 'Materials/cond_inner/prop_values Materials/cond_outer/prop_values
Postprocessors/heat_source/scale_factor
Materials/thermal_strain_inner/thermal_expansion_coeff Materials/thermal_strain_outer/thermal_expansion_coeff
Materials/elasticity_tensor_inner/youngs_modulus Materials/elasticity_tensor_outer/youngs_modulus
Materials/elasticity_tensor_inner/poissons_ratio Materials/elasticity_tensor_outer/poissons_ratio'
to_control = 'stochastic'
check_multiapp_execute_on = false
[]
[data]
type = SamplerPostprocessorTransfer
multi_app = sub
sampler = sample
to_vector_postprocessor = storage
from_postprocessor = 'temp_center_inner  temp_center_outer  temp_end_inner  temp_end_outer
dispx_center_inner dispx_center_outer dispx_end_inner dispx_end_outer
dispz_inner dispz_outer'
[]
[]

[VectorPostprocessors]
[storage]
type = StochasticResults
[]
[]

[Trainers]
[temp_center_inner]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 4
sampler = sample
results_vpp = storage
results_vector = 'data:temp_center_inner'
[]
[temp_center_outer]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 4
sampler = sample
results_vpp = storage
results_vector = 'data:temp_center_outer'
[]
[temp_end_inner]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 4
sampler = sample
results_vpp = storage
results_vector = 'data:temp_end_inner'
[]
[temp_end_outer]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 4
sampler = sample
results_vpp = storage
results_vector = 'data:temp_end_outer'
[]
[dispx_center_inner]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 4
sampler = sample
results_vpp = storage
results_vector = 'data:dispx_center_inner'
[]
[dispx_center_outer]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 4
sampler = sample
results_vpp = storage
results_vector = 'data:dispx_center_outer'
[]
[dispx_end_inner]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 4
sampler = sample
results_vpp = storage
results_vector = 'data:dispx_end_inner'
[]
[dispx_end_outer]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 4
sampler = sample
results_vpp = storage
results_vector = 'data:dispx_end_outer'
[]
[dispz_inner]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 4
sampler = sample
results_vpp = storage
results_vector = 'data:dispz_inner'
[]
[dispz_outer]
type = PolynomialChaosTrainer
execute_on = timestep_end
order = 4
sampler = sample
results_vpp = storage
results_vector = 'data:dispz_outer'
[]
[]

[Outputs]
[out]
type = SurrogateTrainerOutput
trainers = 'temp_center_inner  temp_center_outer  temp_end_inner  temp_end_outer
dispx_center_inner dispx_center_outer dispx_end_inner dispx_end_outer
dispz_inner dispz_outer'
execute_on = FINAL
[]
[]

(modules/combined/examples/stochastic/poly_chaos_train_uniform.i)

Listing 4: Polynomial chaos evaluation input file

[StochasticTools]
[]

[Distributions]
[cond_inner]
type = Uniform
lower_bound = 20
upper_bound = 30
[]
[cond_outer]
type = Uniform
lower_bound = 90
upper_bound = 110
[]
[heat_source]
type = Uniform
lower_bound = 9000
upper_bound = 11000
[]
[alpha_inner]
type = Uniform
lower_bound = 1e-6
upper_bound = 3e-6
[]
[alpha_outer]
type = Uniform
lower_bound = 5e-7
upper_bound = 1.5e-6
[]
[ymod_inner]
type = Uniform
lower_bound = 2e5
upper_bound = 2.2e5
[]
[ymod_outer]
type = Uniform
lower_bound = 3e5
upper_bound = 3.2e5
[]
[prat_inner]
type = Uniform
lower_bound = 0.29
upper_bound = 0.31
[]
[prat_outer]
type = Uniform
lower_bound = 0.19
upper_bound = 0.21
[]
[]

[Samplers]
[sample]
type = MonteCarlo
num_rows = 100000
distributions = 'cond_inner cond_outer heat_source alpha_inner alpha_outer ymod_inner ymod_outer prat_inner prat_outer'
execute_on = INITIAL
[]
[]

[Surrogates]
[temp_center_inner]
type = PolynomialChaos
filename = 'poly_chaos_train_uniform_out_temp_center_inner.rd'
[]
[temp_center_outer]
type = PolynomialChaos
filename = 'poly_chaos_train_uniform_out_temp_center_outer.rd'
[]
[temp_end_inner]
type = PolynomialChaos
filename = 'poly_chaos_train_uniform_out_temp_end_inner.rd'
[]
[temp_end_outer]
type = PolynomialChaos
filename = 'poly_chaos_train_uniform_out_temp_end_outer.rd'
[]
[dispx_center_inner]
type = PolynomialChaos
filename = 'poly_chaos_train_uniform_out_dispx_center_inner.rd'
[]
[dispx_center_outer]
type = PolynomialChaos
filename = 'poly_chaos_train_uniform_out_dispx_center_outer.rd'
[]
[dispx_end_inner]
type = PolynomialChaos
filename = 'poly_chaos_train_uniform_out_dispx_end_inner.rd'
[]
[dispx_end_outer]
type = PolynomialChaos
filename = 'poly_chaos_train_uniform_out_dispx_end_outer.rd'
[]
[dispz_inner]
type = PolynomialChaos
filename = 'poly_chaos_train_uniform_out_dispz_inner.rd'
[]
[dispz_outer]
type = PolynomialChaos
filename = 'poly_chaos_train_uniform_out_dispz_outer.rd'
[]
[]

[VectorPostprocessors]
[stats]
type = PolynomialChaosStatistics
pc_name = 'temp_center_inner  temp_center_outer  temp_end_inner  temp_end_outer
dispx_center_inner dispx_center_outer dispx_end_inner dispx_end_outer
dispz_inner dispz_outer'
compute = 'mean stddev'
[]
[sobol]
type = PolynomialChaosSobolStatistics
pc_name = 'temp_center_inner  temp_center_outer  temp_end_inner  temp_end_outer
dispx_center_inner dispx_center_outer dispx_end_inner dispx_end_outer
dispz_inner dispz_outer'
sensitivity_order = 'first second total'
[]
[storage]
type = EvaluateSurrogate
sampler = sample
model = 'temp_center_inner  temp_center_outer  temp_end_inner  temp_end_outer
dispx_center_inner dispx_center_outer dispx_end_inner dispx_end_outer
dispz_inner dispz_outer'
[]
[]

[Outputs]
csv = true
execute_on = TIMESTEP_END
[]

(modules/combined/examples/stochastic/poly_chaos_uniform.i)

### Statistics

Table 4 shows the statistical results of sampling the thermomechanis model and the polynomial chaos surrogate. and represent the mean and standard deviation of the QoI, and CI is the confidence interval. Note that the confidence interval for the polynomial chaos statistics is not relevant since these values were found analytically using integration techniques. Figure 9 to Figure 11 compares several of the probability distributions of the QoIs between sampling the full-order model and the polynomial chaos surrogate.

Table 4: Statistics Results

QoI95% CI95% CIPC – PC –
609.95(609.85,
610.04)
18.23(18.19,
18.28)
609.9718.23
586.84(586.75,
586.93)
16.63(16.59,
16.67)
586.8516.64
506.31(506.25,
506.38)
12.04(12.01,
12.06)
506.3112.05
507.92(507.86,
507.99)
12.11(12.08,
12.14)
507.9212.12
4.034E-04(4.030E-04,
4.039E-04)
8.626E-05(8.599E-05,
8.653E-05)
4.032E-048.625E-05
4.999E-04(4.993E-04,
5.004E-04)
1.080E-04(1.077E-04,
1.084E-04)
4.996E-041.081E-04
4.223E-04(4.217E-04,
4.230E-04)
1.254E-04(1.251E-04,
1.257E-04)
4.220E-041.256E-04
4.797E-04(4.790E-04,
4.804E-04)
1.352E-04(1.349E-04,
1.356E-04)
4.794E-041.354E-04
6.143E-04(6.135E-04,
6.150E-04)
1.467E-04(1.463E-04,
1.471E-04)
6.139E-041.469E-04
4.997E-04(4.992E-04,
5.003E-04)
1.071E-04(1.067E-04,
1.074E-04)
4.995E-041.071E-04

Figure 9: Histogram

Figure 10: Histogram

Figure 11: Histogram

### Sobol Sensitivities

Sobol sensitivities, or Sobol indicies, are a metric to compare the global sensitivity a parameter has on a QoI. This examples demonstrates several different types of the sensitivities. The first is total sensitivity, which measure the total sensitivity from a parameter, Figure 3 shows these values for each QoI and parameter. The second is a correlative sensitivity, which measures the sensitivity due to a combination of parameters, Figure 12 to Figure 14 show heat maps of these values for several QoIs.

Figure 3: Total Sobol sensitivities

Figure 12: Sobol sensitivities

Figure 13: Sobol sensitivities

Figure 14: Sobol sensitivities