Parameter Study on a Highly Nonlinear Problem
This example assumes that the reader has already visited the example in Parameter Study and is familiar with the fundamental blocks used in parent input files. In this example, the effect of varying the distribution of the uncertain parameters on the distribution of the Quantities of Interest (QoIs) is showcased as well.
Problem Description
The strong formulation of the problem is taken from Chaturantabut and Sorensen (2010) and can be written as
(1)where is a scalar field variable, are the physical coordinates, while and are uncertain parameters with known (or assumed) probability distributions. This equation is supplemented with homogeneous Dirichlet boundary conditions on every side of the square.
Solution of the Problem
To be able to perform a parameter study, the application has to be able to solve the problem with fixed parameters first. The input file used for this purpose is provided in Listing 1. The nominal values of the uncertain parameters are and =9 in this case. There are two blocks in the input file that are worth examining in detail. The first is the Kernels block that shows that a custom test kernel has been implemented to be able to handle the exponential reaction term in Eq. (1). To use this kernel, the user has to add an additional argument for the Stochastic Tools executioner as follows:
cd moose/modules/stochastic_tools/examples/parameter_study/nonlin_diff_react
../../../stochastic_tools-opt -i nonlin_diff_react_sub.i --allow-test-objects
The second atypical block is Controls which is necessary to set up a channel to receive and substitute new parameter samples from the parent application. As shown in the Postprocessors block, the Quantities of Interest (QoIs) are the maximum value (), minimum value () and the average value () of the scalar field variable .
Listing 1: Complete input file for the nonlinear problem using the nominal values of the uncertain parameters.
[Functions<<<{"href": "../../../syntax/Functions/index.html"}>>>]
[source]
type = ParsedFunction<<<{"description": "Function created by parsing a string", "href": "../../../source/functions/MooseParsedFunction.html"}>>>
expression<<<{"description": "The user defined function."}>>> = "100 * sin(2 * pi * x) * sin(2 * pi * y)"
[]
[]
[Mesh<<<{"href": "../../../syntax/Mesh/index.html"}>>>]
[gen]
type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 2
nx<<<{"description": "Number of elements in the X direction"}>>> = 50
xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = 0
xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 1
ny<<<{"description": "Number of elements in the Y direction"}>>> = 50
ymin<<<{"description": "Lower Y Coordinate of the generated mesh"}>>> = 0
ymax<<<{"description": "Upper Y Coordinate of the generated mesh"}>>> = 1
[]
[]
[Variables<<<{"href": "../../../syntax/Variables/index.html"}>>>]
[U]
family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = lagrange
order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = first
[]
[]
[Kernels<<<{"href": "../../../syntax/Kernels/index.html"}>>>]
[diffusion]
type = Diffusion<<<{"description": "The Laplacian operator ($-\\nabla \\cdot \\nabla u$), with the weak form of $(\\nabla \\phi_i, \\nabla u_h)$.", "href": "../../../source/kernels/Diffusion.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = U
[]
[nonlin_function]
type = ExponentialReaction
variable = U
mu1 = 0.3
mu2 = 9
[]
[source]
type = BodyForce<<<{"description": "Demonstrates the multiple ways that scalar values can be introduced into kernels, e.g. (controllable) constants, functions, and postprocessors. Implements the weak form $(\\psi_i, -f)$.", "href": "../../../source/kernels/BodyForce.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = U
function<<<{"description": "A function that describes the body force"}>>> = source
[]
[]
[BCs<<<{"href": "../../../syntax/BCs/index.html"}>>>]
[dirichlet_all]
type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../source/bcs/DirichletBC.html"}>>>
variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = U
boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'left right top bottom'
value<<<{"description": "Value of the BC"}>>> = 0
[]
[]
[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
type = Steady
solve_type = PJFNK
petsc_options_iname = '-pc_type -pc_hypre_type'
petsc_options_value = 'hypre boomeramg'
[]
[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
[max]
type = ElementExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../../../source/postprocessors/ElementExtremeValue.html"}>>>
variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = U
[]
[min]
type = ElementExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../../../source/postprocessors/ElementExtremeValue.html"}>>>
variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = U
value_type<<<{"description": "Type of extreme value to return. 'max' returns the maximum value. 'min' returns the minimum value. 'max_abs' returns the maximum of the absolute value."}>>> = min
[]
[average]
type = ElementAverageValue<<<{"description": "Computes the volumetric average of a variable", "href": "../../../source/postprocessors/ElementAverageValue.html"}>>>
variable<<<{"description": "The name of the variable that this object operates on"}>>> = U
[]
[]
[Controls<<<{"href": "../../../syntax/Controls/index.html"}>>>]
[stochastic]
type = SamplerReceiver<<<{"description": "Control for receiving data from a Sampler via SamplerParameterTransfer.", "href": "../../../source/controls/SamplerReceiver.html"}>>>
[]
[]
[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
[]Parent application Input
As described in Parameter Study in detail, one needs a driver input (or parent input) to perform a parameter study. Two parent input files are provided for this example in Listing 2 and Listing 3. The first considers the uncertain parameters to be uniformly distributed around their nominal values , , while the second one assumes normal distribution . The only difference between the two input files is the Distributions block where the assumed probability distributions are defined for the uncertain parameters.
Listing 2: Complete input file for the driver of the parameter study with uniformly distributed uncertain parameters.
[StochasticTools<<<{"href": "../../../syntax/StochasticTools/index.html"}>>>]
[]
[Distributions<<<{"href": "../../../syntax/Distributions/index.html"}>>>]
[mu1]
type = Uniform<<<{"description": "Continuous uniform distribution.", "href": "../../../source/distributions/Uniform.html"}>>>
lower_bound<<<{"description": "Distribution lower bound"}>>> = 0.21
upper_bound<<<{"description": "Distribution upper bound"}>>> = 0.39
[]
[mu2]
type = Uniform<<<{"description": "Continuous uniform distribution.", "href": "../../../source/distributions/Uniform.html"}>>>
lower_bound<<<{"description": "Distribution lower bound"}>>> = 6.3
upper_bound<<<{"description": "Distribution upper bound"}>>> = 11.7
[]
[]
[Samplers<<<{"href": "../../../syntax/Samplers/index.html"}>>>]
[hypercube]
type = LatinHypercube<<<{"description": "Latin Hypercube Sampler.", "href": "../../../source/samplers/LatinHypercubeSampler.html"}>>>
num_rows<<<{"description": "The size of the square matrix to generate."}>>> = 5000
distributions<<<{"description": "The distribution names to be sampled, the number of distributions provided defines the number of columns per matrix."}>>> = 'mu1 mu2'
[]
[]
[MultiApps<<<{"href": "../../../syntax/MultiApps/index.html"}>>>]
[runner]
type = SamplerFullSolveMultiApp<<<{"description": "Creates a full-solve type sub-application for each row of each Sampler matrix.", "href": "../../../source/multiapps/SamplerFullSolveMultiApp.html"}>>>
sampler<<<{"description": "The Sampler object to utilize for creating MultiApps."}>>> = hypercube
input_files<<<{"description": "The input file for each App. If this parameter only contains one input file it will be used for all of the Apps. When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = 'nonlin_diff_react_sub.i'
mode<<<{"description": "The operation mode, 'normal' creates one sub-application for each row in the Sampler and 'batch-reset' and 'batch-restore' creates N sub-applications, where N is the minimum of 'num_rows' in the Sampler and floor(number of processes / min_procs_per_app). To run the rows in the Sampler, 'batch-reset' will destroy and re-create sub-apps as needed, whereas the 'batch-restore' will backup and restore sub-apps to the initial state prior to execution, without destruction."}>>> = batch-restore
[]
[]
[Transfers<<<{"href": "../../../syntax/Transfers/index.html"}>>>]
[parameters]
type = SamplerParameterTransfer<<<{"description": "Copies Sampler data to a SamplerReceiver object.", "href": "../../../source/transfers/SamplerParameterTransfer.html"}>>>
to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = runner
sampler<<<{"description": "A the Sampler object that Transfer is associated.."}>>> = hypercube
parameters<<<{"description": "A list of parameters (on the sub application) to control with the Sampler data. The order of the parameters listed here should match the order of the items in the Sampler."}>>> = 'Kernels/nonlin_function/mu1 Kernels/nonlin_function/mu2'
[]
[results]
type = SamplerPostprocessorTransfer<<<{"description": "Transfers data from Postprocessors on the sub-application to a VectorPostprocessor on the master application.", "href": "../../../source/transfers/SamplerPostprocessorTransfer.html"}>>>
from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = runner
sampler<<<{"description": "A the Sampler object that Transfer is associated.."}>>> = hypercube
to_vector_postprocessor<<<{"description": "The name of the VectorPostprocessor in the MultiApp to transfer values to."}>>> = results
from_postprocessor<<<{"description": "The name(s) of the Postprocessor(s) on the sub-app to transfer from."}>>> = 'max min average'
[]
[]
[VectorPostprocessors<<<{"href": "../../../syntax/VectorPostprocessors/index.html"}>>>]
[results]
type = StochasticResults<<<{"description": "Storage container for stochastic simulation results coming from a Postprocessor.", "href": "../../../source/vectorpostprocessors/StochasticResults.html"}>>>
[]
[]
[Reporters<<<{"href": "../../../syntax/Reporters/index.html"}>>>]
[stats]
type = StatisticsReporter<<<{"description": "Compute statistical values of a given VectorPostprocessor objects and vectors.", "href": "../../../source/reporters/StatisticsReporter.html"}>>>
vectorpostprocessors<<<{"description": "List of VectorPostprocessor(s) to utilized for statistic computations."}>>> = results
compute<<<{"description": "The statistic(s) to compute for each of the supplied vector postprocessors."}>>> = 'mean'
ci_method<<<{"description": "The method to use for computing confidence level intervals."}>>> = 'percentile'
ci_levels<<<{"description": "A vector of confidence levels to consider, values must be in (0, 1)."}>>> = '0.05'
[]
[]
[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'FINAL'
[]Listing 3: Complete input file for the driver of the parameter study with normally distributed uncertain parameters.
[Distributions<<<{"href": "../../../syntax/Distributions/index.html"}>>>]
[mu1]
type = Normal<<<{"description": "Normal distribution", "href": "../../../source/distributions/Normal.html"}>>>
mean<<<{"description": "Mean (or expectation) of the distribution."}>>> = 0.3
standard_deviation<<<{"description": "Standard deviation of the distribution "}>>> = 0.045
[]
[mu2]
type = Normal<<<{"description": "Normal distribution", "href": "../../../source/distributions/Normal.html"}>>>
mean<<<{"description": "Mean (or expectation) of the distribution."}>>> = 9
standard_deviation<<<{"description": "Standard deviation of the distribution "}>>> = 1.35
[]
[]For the sampling of the uncertain parameters, a Latin Hypercube Sampling (LHS) strategy is utilized. Altogether 5000 parameter samples are created for the model. Furthermore, the parent application is executed in a "batch-restore" mode, which provides a memory-efficient way to run sub-applications. For the comparison between different running modes the interested reader is referred to Stochastic Tools Batch Mode.
The objects in the Transfers block are responsible for the communication between the parent and sub-applications. It streams parameter samples to sub-applications and receives the corresponding values for the selected QoIs. It is visible that in this example the type of the parameter transfer object is SamplerParameterTransfer which streams the parameter samples to a SamplerReceiver object (in Controls block) in the sub-application. This object then plugs the new parameter values into kernels, materials or boundary conditions. Unfortunately, this requires the parameters to be controllable in the sub-application, which might not be true in every case. For this specific example, the controllability of the parameters in ExponentialReaction kernel is ensured by the last two commands in the validParams function.
InputParameters
ExponentialReaction::validParams()
{
InputParameters params = Kernel::validParams();
params.addClassDescription(
"Implements a simple reaction term with the following weak form "
"$(\\psi_i,\\frac{\\mu_1}{\\mu_2}\\left[e^{\\mu_2 u_h}-1\\right])$.");
params.addParam<Real>("mu1", 0.3, "First coefficient in the nonlinear term.");
params.addParam<Real>("mu2", 9, "Second coefficient in the nonlinear term.");
params.declareControllable("mu1");
params.declareControllable("mu2");
return params;
}
If the target parameters are not controllable, one can use a command line based communication between parent and sub-applications. For more information about this approach see the example covered in Polynomial Chaos Surrogate.
To run the parent application, it is still necessary to enable test objects using the following command
cd moose/modules/stochastic_tools/examples/parameter_study/nonlin_diff_react
../../../stochastic_tools-opt -i nonlin_diff_react_parent_uniform.i --allow-test-objects
Stochastic Results
The distributions of the QoIs for the uniformly distributed uncertain parameters are presented in Figure 1, Figure 2 and Figure 3. The same distributions with normally distributed uncertain parameters are shown in Figure 4, Figure 5 and Figure 6.
The estimated mean values of the QoIs with the corresponding confidence intervals are presented below assuming uniformly distributed parameters:
The same statistics for the normally distributed parameters are the following:
Figure 1: Resulting distribution of with uniformly distributed parameters.
Figure 2: Resulting distribution of with uniformly distributed parameters.
Figure 3: Resulting distribution of with uniformly distributed parameters.
Figure 4: Resulting distribution of with normally distributed parameters.
Figure 5: Resulting distribution of with normally distributed parameters.
Figure 6: Resulting distribution of with normally distributed parameters.
References
- Saifon Chaturantabut and Danny C Sorensen.
Nonlinear model reduction via discrete empirical interpolation.
SIAM Journal on Scientific Computing, 32(5):2737–2764, 2010.[BibTeX]