POD Reduced Basis Surrogate

This example is meant to demonstrate how a POD reduced basis surrogate model is trained and used on a parametric problem.

Problem statement

The full-order model is a one energy group, fixed-source diffusion-reaction problem, adopted from Prince and Ragusa (2019). The geometry for this problem is presented in Figure 1. The problem has four different material regions, from which three (1, 2 and 3) act as fixed sources.

Figure 1: The geometry of the diffusion-reaction problem used in this example (Prince and Ragusa (2019)).

The fixed-source diffusion-reaction problem with space dependent coefficients can be expressed as:

−∇⋅[D(r)∇ψ(r)]+Σa(r)ψ(r)=q(r) ,r∈Ω-\nabla \cdot\left[D(\textbf{r})\nabla\psi(\textbf{r})\right]+\Sigma_a(\textbf{r})\psi(\textbf{r}) = q(\textbf{r}) \,, \quad \textbf{r}\in\Omega(1)

where D(r)D(\textbf{r}) is the diffusion coefficient, Σa(r)\Sigma_a(\textbf{r}) is the reaction coefficient, q(r)q(\textbf{r}) is the fixed source term and field variable ψ(r)\psi(\textbf{r}) is the solution of interest. Furthermore, Ω\Omega denotes the internal domain, without the boundaries, which can be partitioned into four sub-domains corresponding to the four material regions (Ω1\Omega_1, Ω2\Omega_2, Ω3\Omega_3 and Ω4\Omega_4). This equation needs to be supplemented with boundary conditions. For the symmetry lines (dashed lines in Figure 1, denoted by ∂Ωsym\partial\Omega_{sym}) a homogeneous Neumann condition is used:

−D(r)∇ψ(r)⋅n(r)=0 ,r∈∂Ωsym-D(\textbf{r})\nabla\psi(\textbf{r})\cdot \textbf{n}(\textbf{r}) = 0 \,, \quad \textbf{r}\in\partial\Omega_{sym}

while the rest of the boundaries (in the reflector, denoted by ∂Ωrefl\partial\Omega_{refl}) are treated with homogeneous Dirichlet conditions:

ψ(r)=0 ,r∈∂Ωrefl\psi(\textbf{r}) = 0 \,, \quad \textbf{r}\in\partial\Omega_{refl}

This problem is parametric in a sense that the solution depends on the values of the coefficients and the source: ψ=ψ(r,D(r),Σa(r),q(r))\psi = \psi\left(\textbf{r},D(\textbf{r}),\Sigma_a(\textbf{r}),q(\textbf{r})\right). In this example, material region-wise constant coefficients and source terms are considered yielding eight uncertain parameters altogether (assuming that Region 4 does not have a source). The material properties in each region have Uniform distributions (U(a,b)\mathcal{U}(a,b)) specified in Table 1 with aa and bb being the lower and upper bounds.

Table 1: The distributions of the uncertain parameters used in this problem (Prince and Ragusa (2019)).

ParameterSymbolDistribution
Diffusion coefficient in Region 1 (cm)\left(cm\right)D1D_1∼U(0.2,0.8)\sim\mathcal{U}(0.2, 0.8)
Diffusion coefficient in Region 2 (cm)\left(cm\right)D2D_2∼U(0.2,0.8)\sim\mathcal{U}(0.2, 0.8)
Diffusion coefficient in Region 3 (cm)\left(cm\right)D3D_3∼U(0.2,0.8)\sim\mathcal{U}(0.2, 0.8)
Diffusion coefficient in Region 4 (cm)\left(cm\right)D4D_4∼U(0.15,0.6)\sim\mathcal{U}(0.15, 0.6)
Reaction coefficient in Region 1 (cm−1)\left(cm^{-1}\right)Σa,1\Sigma_{a,1}∼U(0.0425,0.17)\sim\mathcal{U}(0.0425, 0.17)
Reaction coefficient in Region 2 (cm−1)\left(cm^{-1}\right)Σa,2\Sigma_{a,2}∼U(0.065,0.26)\sim\mathcal{U}(0.065, 0.26)
Reaction coefficient in Region 3 (cm−1)\left(cm^{-1}\right)Σa,3\Sigma_{a,3}∼U(0.04,0.16)\sim\mathcal{U}(0.04, 0.16)
Reaction coefficient in Region 4 (cm−1)\left(cm^{-1}\right)Σa,4\Sigma_{a,4}∼U(0.005,0.02)\sim\mathcal{U}(0.005, 0.02)
Fixed-source in Region 1 (ncm3⋅s)\left(\frac{n}{cm^3 \cdot s}\right)q1q_1∼U(5,20)\sim\mathcal{U}(5, 20)
Fixed-source in Region 2 (ncm3⋅s)\left(\frac{n}{cm^3 \cdot s}\right)q2q_2∼U(5,20)\sim\mathcal{U}(5, 20)
Fixed-source in Region 3 (ncm3⋅s)\left(\frac{n}{cm^3 \cdot s}\right)q3q_3∼U(5,20)\sim\mathcal{U}(5, 20)

It is important to mention that POD-RB surrogates are only efficient when the original problem has an affine decomposition. For more information about affine decomposition see PODReducedBasisTrainer. Luckily, the problem at hand has an affine decomposition in the following form:

−∑i=14∇⋅[Di(r)∇ψ(r)]+∑i=14Σa,i(r)ψ(r)=∑i=13qi(r) ,r∈Ω-\sum\limits_{i=1}^{4}\nabla \cdot\left[D_i(\textbf{r})\nabla\psi(\textbf{r})\right]+\sum\limits_{i=1}^{4}\Sigma_{a,i}(\textbf{r})\psi(\textbf{r}) = \sum\limits_{i=1}^{3}q_i(\textbf{r}) \,, \quad \textbf{r}\in\Omega

where Di(r)D_i(\textbf{r}), Σa,i(r)\Sigma_{a,i}(\textbf{r}) and qi(r)q_i(\textbf{r}) take the values of DiD_i, Σa,i\Sigma_{a,i} and qiq_i when r∈Ωi\textbf{r}\in\Omega_i and 0 otherwise.

Solving the problem without uncertain parameters

The first step towards creating a POD-RB surrogate model is the generation of a full-order problem which can solve Eq. (1) with fixed parameters. There are three important factors that need to be considered while preparing the input file for this problem:

  1. The user must specify vector tags in the Problem block for each component in the affine decomposition of the system. In this example, as shown in Listing 1, 8 vector tags are specified for the eight uncertain parameters. These do not introduce extra work in the full-order model, however they help to identify the affine components throughout the training phase.

    Listing 1: Complete input file for the heat equation problem in this study.

    [Problem]
  type = FEProblem
  extra_tag_vectors = 'diff0 diff1 diff2 diff3 abs0 abs1 abs2 abs3 src0 src1 src2'
[]
    (contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/sub.i)

  2. The input file has to reflect the affine decomposition of the problem. This means that the Kernels, BCs and Materials have to be created in a way that they correspond to the components in the affine decomposition. This is shown in Listing 2. Note that a separate kernel has been created for every single term in the decomposition. The vector tags in the Problem block are then applied to these kernels to ensure that the affine components are correctly identified throughout the training phase.

    Listing 2: Complete input file for the heat equation problem in this study.

    [Kernels]
  [diff0]
    type = MatDiffusion
    variable = psi
    diffusivity = D0
    extra_vector_tags = 'diff0'
    block = 0
  []
  [diff1]
    type = MatDiffusion
    variable = psi
    diffusivity = D1
    extra_vector_tags = 'diff1'
    block = 1
  []
  [diff2]
    type = MatDiffusion
    variable = psi
    diffusivity = D2
    extra_vector_tags = 'diff2'
    block = 2
  []
  [diff3]
    type = MatDiffusion
    variable = psi
    diffusivity = D3
    extra_vector_tags = 'diff3'
    block = 3
  []
  [abs0]
    type = MaterialReaction
    variable = psi
    coefficient = absxs0
    extra_vector_tags = 'abs0'
    block = 0
  []
  [abs1]
    type = MaterialReaction
    variable = psi
    coefficient = absxs1
    extra_vector_tags = 'abs1'
    block = 1
  []
  [abs2]
    type = MaterialReaction
    variable = psi
    coefficient = absxs2
    extra_vector_tags = 'abs2'
    block = 2
  []
  [abs3]
    type = MaterialReaction
    variable = psi
    coefficient = absxs3
    extra_vector_tags = 'abs3'
    block = 3
  []
  [src0]
    type = BodyForce
    variable = psi
    value = 1
    extra_vector_tags = 'src0'
    block = 0
  []
  [src1]
    type = BodyForce
    variable = psi
    value = 1
    extra_vector_tags = 'src1'
    block = 1
  []
  [src2]
    type = BodyForce
    variable = psi
    value = 1
    extra_vector_tags = 'src2'
    block = 2
  []
[]

[Materials]
  [D0]
    type = GenericConstantMaterial
    prop_names = D0
    prop_values = 1
    block = 0
  []
  [D1]
    type = GenericConstantMaterial
    prop_names = D1
    prop_values = 1
    block = 1
  []
  [D2]
    type = GenericConstantMaterial
    prop_names = D2
    prop_values = 1
    block = 2
  []
  [D3]
    type = GenericConstantMaterial
    prop_names = D3
    prop_values = 1
    block = 3
  []
  [absxs0]
    type = GenericConstantMaterial
    prop_names = absxs0
    prop_values = 1
    block = 0
  []
  [absxs1]
    type = GenericConstantMaterial
    prop_names = absxs1
    prop_values = 1
    block = 1
  []
  [absxs2]
    type = GenericConstantMaterial
    prop_names = absxs2
    prop_values = 1
    block = 2
  []
  [absxs3]
    type = GenericConstantMaterial
    prop_names = absxs3
    prop_values = 1
    block = 3
  []
[]
    (contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/sub.i)

  3. The values for each uncertain parameter have to be set to 1 by default. This is necessary because the PODReducedBasisTrainer uses the same input file to create the affine constituent operators. This ensures that the mentioned operators are not influenced by the parameter-dependent multipliers. Of course, these values are not fixed and are changed by the main application throughout the simulations to values aligned with those in Table 1. However, the default values in the input file should be set to one. This is shown in Listing 2 as well.

Training surrogate models

For the details of the training procedure of a POD-RB surrogate, see PODReducedBasisTrainer. The first step is the collection of data, which involves the repeated execution of the full-order problem with different parameter combinations and the saving of the full solution vectors. These solution vectors are often referred to as snapshots and this naming is preferred in this example as well. This step is managed by the main input file which creates parameter samples, transfers them to the sub-application and collects the results from the completed computations.

The snapshot collection phase starts with the definition of the distributions in the Distributions block. The uniform distributions for all the 8 uncertain parameters are specified as:

[Distributions]
  [D012_dist]
    type = Uniform
    lower_bound = 0.2
    upper_bound = 0.8
  []
  [D1_dist]
    type = Uniform
    lower_bound = 0.2
    upper_bound = 0.8
  []
  [D2_dist]
    type = Uniform
    lower_bound = 0.2
    upper_bound = 0.8
  []
  [D3_dist]
    type = Uniform
    lower_bound = 0.15
    upper_bound = 0.6
  []
  [absxs0_dist]
    type = Uniform
    lower_bound = 0.0425
    upper_bound = 0.17
  []
  [absxs1_dist]
    type = Uniform
    lower_bound = 0.065
    upper_bound = 0.26
  []
  [absxs2_dist]
    type = Uniform
    lower_bound = 0.04
    upper_bound = 0.16
  []
  [absxs3_dist]
    type = Uniform
    lower_bound = 0.005
    upper_bound = 0.02
  []
  [src_dist]
    type = Uniform
    lower_bound = 5
    upper_bound = 20
  []
[]
(contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/trainer.i)

As a next step, the underlying distributions are sampled to create parameter combinations. This is done using a LatinHypercube defined in the Samplers block. It is visible that 100 samples are prepared, meaning that 100 snapshots will be collected for the generation of the surrogates.

[Samplers]
  [sample]
    type = LatinHypercube
    distributions = 'D012_dist D012_dist D012_dist D3_dist
                     absxs0_dist absxs1_dist absxs2_dist absxs3_dist
                     src_dist src_dist src_dist'
    num_rows = 100
    execute_on = PRE_MULTIAPP_SETUP
    max_procs_per_row = 1
  []
[]
(contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/trainer.i)

To be able to create the reduced operators for the surrogate model, a custom MultiApp, PODFullSolveMultiApp, has been created. This object is responsible for executing sub-problems using different combinations of parameter values provided by the sampler. The secondary function of this object is to create the action of the full-order operators on the basis functions of the reduced subspace. Therefore, this object has to be executed twice in the same simulation. It is visible that unlike a regular SamplerFullSolveMultiApp, this custom object has to know certain parameters of the trainer as well.

[MultiApps]
  [sub]
    type = PODFullSolveMultiApp
    input_files = sub.i
    sampler = sample
    trainer_name = 'pod_rb'
    execute_on = 'timestep_begin final'
    max_procs_per_app = 1
  []
[]
(contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/trainer.i)

In terms of the Transfers block, besides sending the actual parameter samples to the sub-applications, in this intrusive procedure, the snapshots need to be collected from the sub-applications, the basis functions need to be sent back to different sub-applications and the action of the operators on the basis functions need to be collected as well. This requires four transfer objects. The two custom types (PODSamplerSolutionTransfer and PODResidualTransfer) are specifically used to support PODReducedBasisTrainer at this moment.

[Transfers]
  [param]
    type = SamplerParameterTransfer
    to_multi_app = sub
    sampler = sample
    parameters = 'Materials/D0/prop_values
                  Materials/D1/prop_values
                  Materials/D2/prop_values
                  Materials/D3/prop_values
                  Materials/absxs0/prop_values
                  Materials/absxs1/prop_values
                  Materials/absxs2/prop_values
                  Materials/absxs3/prop_values
                  Kernels/src0/value
                  Kernels/src1/value
                  Kernels/src2/value'
    execute_on = 'timestep_begin'
    check_multiapp_execute_on = false
  []
  [data]
    type = PODSamplerSolutionTransfer
    from_multi_app = sub
    sampler = sample
    trainer_name = 'pod_rb'
    execute_on = 'timestep_begin'
    check_multiapp_execute_on = false
  []
  [mode]
    type = PODSamplerSolutionTransfer
    to_multi_app = sub
    sampler = sample
    trainer_name = 'pod_rb'
    execute_on = 'final'
    check_multiapp_execute_on = false
  []
  [res]
    type = PODResidualTransfer
    from_multi_app = sub
    sampler = sample
    trainer_name = 'pod_rb'
    execute_on = 'final'
    check_multiapp_execute_on = false
  []
[]
(contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/trainer.i)

Next, the PODReducedBasisTrainer is set up in the Trainers block. The trainer stores the snapshots and uses them to create basis functions for the reduced subspaces. Furthermore, it is also responsible for creating the reduced-order operators, therefore it needs to be executed twice in the training process. The trainer object needs to know what variable needs to be reduced and the names of the vector tags from the sub-application to be able to identify the affine constituent operators. Furthermore, using the "tag_types" input argument, the user has to specify if the reduced affine constituent operator acts on the variable or not. The ordering must be the same as the names of the vector tags. The meaning of the energy retention limits is discussed in PODReducedBasisTrainer.

[Trainers]
  [pod_rb]
    type = PODReducedBasisTrainer
    var_names = 'psi'
    error_res = '1e-9'
    tag_names = 'diff0 diff1 diff2 diff3 abs0 abs1 abs2 abs3 src0 src1 src2'
    tag_types = 'op op op op op op op op src src src'
    execute_on = 'timestep_begin final'
  []
[]
(contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/trainer.i)

As a last step in the training process, the basis functions, reduced operators and every necessary information for the surrogate are saved into an .rd file. This file can be then used to construct a surrogate model without the need to repeat the training process.

[Outputs]
  [out]
    type = SurrogateTrainerOutput
    trainers = 'pod_rb'
    execute_on = FINAL
  []
[]
(contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/trainer.i)

Evaluation of surrogate models

To evaluate surrogate models, a new main input file has to be created. In this example, the same distributions are defined for the parameters as used in the training phase. Therefore, the content of the Distributions block is identical to the one in the trainer input file. As a next step, new samples are generated using these distributions. Again, a LatinHypercube is employed for this task, however this time the number of samples is increased to 1000 since the surrogates are orders of magnitudes faster than the full-order model.

[Samplers]
  [sample]
    type = LatinHypercube
    distributions = 'D012_dist D012_dist D012_dist D3_dist
                     absxs0_dist absxs1_dist absxs2_dist absxs3_dist
                     src_dist src_dist src_dist'
    num_rows = 1000
    execute_on = PRE_MULTIAPP_SETUP
  []
[]
(contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/surr.i)

A PODReducedBasisSurrogate is created in the Surrogates block. It is constructed using the information available within the corresponding .rd file and allows the user to change of the rank of the sub-spaces used for different variables through "change_rank" and "new_ranks" parameters.

[Surrogates]
  [rbpod]
    type = PODReducedBasisSurrogate
    filename = 'trainer_out_pod_rb.rd'
    change_rank = 'psi'
    new_ranks = '40'
  []
[]
(contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/surr.i)

These surrogate models can be evaluated at the points defined in the testing sample batch by a PODSurrogateTester object in the VectorPostprocessors block. In this case, the Quantity of Interest (QoI) is the nodal L2L^2 norm of the solution for ψ\psi.

[VectorPostprocessors]
  [res]
    type = PODSurrogateTester
    model = rbpod
    sampler = sample
    variable_name = 'psi'
    to_compute = nodal_l2
  []
[]
(contrib/moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg/surr.i)

Running the input files

Since the sub-applications use test objects, one has to allow the executioner to use them by specifying the following argument on the command line:

cd moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg
../../../../stochastic_tools-opt -i trainer.i --allow-test-objects

The same is true for the surrogate input file as well, therefore one needs to start the executions as follows:

cd moose/modules/stochastic_tools/examples/surrogates/pod_rb/2d_multireg
../../../../stochastic_tools-opt -i surr.i --allow-test-objects

Results and Analysis

In the following subsection a short analysis is provided for the results obtained using the example input files. As already mentioned, the problem has 8 uncertain parameters and altogether 100 parameter samples are generated using LatinHypercube to obtain snapshots for the training. Three examples of the snapshots are presented in Figure 6, Figure 7 and Figure 8. It is visible that depending on the actual parameter combination, the profile of the solution can change considerably.

Figure 6: Snap. example #1.

Figure 7: Snap. example #2.

Figure 8: Snap. example #3.

After all of the snapshots are obtained, the basis functions of the reduced subspaces are extracted. In this scenario, an energy retention limit of 0.999 999 999 is used in the trainer which will keep 55 basis functions for the reduced subspace. The decay of the eigenvalues of the snapshot correlation matrix is shown in Figure 2. The reduced operators are then computed using these 55 basis functions.

102030405010010k1M100M10B
Basis functionCorresponding eigenvalue

Figure 2: Scree plot of the eigenvalues of the correlation matrix.

As a next step, two surrogate models are prepared using the "change_rank" and "new_ranks" parameters of PODReducedBasisSurrogate to change the size of the reduced system. The first surrogate model has 1 basis function, while the other has 8. Both models are then run on a 1000 sample test set and the nodal L2L^2 norms of the approximate solutions are saved. Additionally, a full-order model was executed on the same test set and the results are saved to serve as reference values. Figure 3 presents the results with the surrogate model built with 1 basis function only. It is visible that the distribution of the QoI (nodal L2L^2 norm) on the test set is considerably different than the reference distribution.

Figure 3: The histogram of the QoI for the full-order reference model and the surrogate built with 1 basis function.

Figure 4 shows the distribution of the QoI obtained by a surrogate with 8 basis functions. It is visible that the difference between the reference values and those from the surrogate is negligible.

Figure 4: The histogram of the QoI for the full-order reference model and the surrogate built with 8 basis function.

To see the convergence of the results from the surrogate to those of the full-order model, the surrogate model is run multiple times with different ranks and the following error indicators are computed for each sample in the test set:

ei=∣∣ψRef−ψSurr∣∣l2∣∣ψRef∣∣l2.i=1,...,1000e_i = \frac{||\psi_{Ref}-\psi_{Surr}||_{l^2}}{||\psi_{Ref}||_{l^2}}.\quad i=1,...,1000

Then, the maximum and average relative errors are recorder as function of the number of basis functions used. Figure 5 shows the results. It is visible that by increasing the number of basis functions, both error indicators decrease rapidly.

9123456789102310μ100μ0.0010.010.11
Average eMaximum eNumber of basis functions usedRelative error

Figure 5: The convergence of averaged quantities of interest.

Lastly, the computation time full-order model on the test set is compared to the combined cost of training and evaluating a POD-RB surrogate model in Table 2. The test has been carried out on one processor only, not using the parallel capabilities of the MultiApp system. The results show that it is beneficial to create a POD-RB surrogate if more than 148 evaluations are needed. This assumes that the full-order evaluation time can be equally distributed among the 1000 test samples (0.779 s/sample). By dividing the training time with this number we get a critical sample number above which the generation of a surrogate model is a better alternative.

Table 2: The computation time of the full-order solutions on the test set compared to the cost of training a surrogate and evaluating it on the same test set.

ProcessExecution time (s)
Evaluation of the full-order model on the 1000 sample test set779.5
Training a POD-RB surrogate using 100 samples116.2
Evaluation of the POD-RB surrogate on the 1000 sample test set (4 basis functions)0.592
Evaluation of the POD-RB surrogate on the 1000 sample test set (8 basis functions)0.937
Evaluation of the POD-RB surrogate on the 1000 sample test set (16 basis functions)1.576

References

  1. Zachary M Prince and Jean C Ragusa. Parametric uncertainty quantification using proper generalized decomposition applied to neutron diffusion. International Journal for Numerical Methods in Engineering, 119(9):899–921, 2019.[BibTeX]