Constrained Optimization Example: Shape Optimization

Problem Description

In this problem, we focus on an annular (ring-shaped) domain subjected to thermal conditions. The domain is defined by two concentric circles with an inner and outer radii. The physical problem is heat conduction governed by the heat equation with convective flux type boundary conditions at the inner radius and an insulated outer radius. The convective heat transfer coefficient is dependent on the inner radius. A constant source term is distributed throughout the domain.

The optimization goal is to find the inner and outer radii that minimize the maximum temperature within the domain while maintaining a fixed volume of the annulus. In other words, the design variables are the inner and outer radii, and the objective function is the maximum temperature within the annulus. This is a constrained optimization problem because the volume of the annulus needs to remain constant during the optimization process.

Example: Shape Optimization

The MOOSE Optimization module provides a flexible framework for solving optimization problems in MOOSE. Shape optimization is a branch of optimization focused on finding the most efficient, effective, or optimal shape for a given system, subject to certain constraints. The goal is usually to minimize or maximize some performance criterion, such as material usage, weight, aerodynamic drag, or structural stress, or as in this case the minimizing the maximum temperature. This page will demonstrate constrained shape optimization of an annulus while satisfying an equality volume constraint.

For constrained optimization using TAO, MOOSE only supports interfacing with the TAOALMM algorithm. The TAOALMM lets user to solve problems in the form below:

$var element = document.getElementById("moose-equation-8c4ccc8c-ba51-42cb-89ac-29acae45dc60");katex.render("\\begin{aligned} \\min_x \\quad & f(x) \\\\ \\textit{s.t.} \\quad & g(x) = 0 \\\\ \\quad & h(x) \\geq 0 \\\\ \\quad & l \\leq x \\leq u \\\\ \\end{aligned}", element, {displayMode:true,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

where

Objective function: $var element = document.getElementById("moose-equation-c3cd0211-15a3-42d3-be3b-9d3cb57343f2");katex.render("f(x)", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$
Equality constraints: $var element = document.getElementById("moose-equation-9809ef61-9f9c-4970-a981-b5169bdceec2");katex.render("g(x) = 0", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$
Inequality constraints: $var element = document.getElementById("moose-equation-15939e5f-5957-4b0e-afda-6a4502d08dfc");katex.render("h(x) \\geq 0", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$
Bounds : $var element = document.getElementById("moose-equation-ca4461db-5ebd-48cc-bfb4-045f861d82a6");katex.render("l \\leq x \\leq u", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$

To solve this problem we use the MultiApps system. For this problem we have a main app that will perform all the optimization and interfacing with TAO, while the forward app will solve the physics and calculate the objective and equality constraints. These values will then be sent using Transfers System to the main app.

Forward Problem Sub-App

In the forward app, a heat conduction problem is solved with convective flux type boundary condition that changes with the inner radius. The outer radius is insulated and there is a constant source term throughout. The objective is find the two radii that will have the lowest max temperature while satisfying a volume constraint.

In the forward app, the two optimization parameters control the inner and outer values are used to displace the boundaries of the ParsedOptimizationFunction to calculate the cartesian boundary displacement. These values are used displace the boundaries of the annulus to the correct position using Dirichlet boundary conditions. A diffusion problem on the displacements then allow for the interior elements to be smoothed. For the temperature field, the physics is solved on the displaced mesh to calculate the current objective.

Listing 1: Forward App Shape Optimization

[Functions<<<{"href": "../../../syntax/Functions/index.html"}>>>]
  [r1_x]
    type = ParsedOptimizationFunction<<<{"description": "Function used for optimization that uses a parsed expression with parameter dependence.", "href": "../../../source/functions/ParsedOptimizationFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = 'r1 * cos((atan(y/x)))'
    param_symbol_names<<<{"description": "Names of parameters in 'expression' being optimized."}>>> = 'r0 r1'
    param_vector_name<<<{"description": "Reporter or VectorPostprocessor vector containing parameter values."}>>> = 'params/radii'
  []
  [r1_y]
    type = ParsedOptimizationFunction<<<{"description": "Function used for optimization that uses a parsed expression with parameter dependence.", "href": "../../../source/functions/ParsedOptimizationFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = 'r1 * sin((atan(y/x)))'
    param_symbol_names<<<{"description": "Names of parameters in 'expression' being optimized."}>>> = 'r0 r1'
    param_vector_name<<<{"description": "Reporter or VectorPostprocessor vector containing parameter values."}>>> = 'params/radii'
  []
  [r0_x]
    type = ParsedOptimizationFunction<<<{"description": "Function used for optimization that uses a parsed expression with parameter dependence.", "href": "../../../source/functions/ParsedOptimizationFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = 'r0 * cos((atan(y/x)))'
    param_symbol_names<<<{"description": "Names of parameters in 'expression' being optimized."}>>> = 'r0 r1'
    param_vector_name<<<{"description": "Reporter or VectorPostprocessor vector containing parameter values."}>>> = 'params/radii'
  []
  [r0_y]
    type = ParsedOptimizationFunction<<<{"description": "Function used for optimization that uses a parsed expression with parameter dependence.", "href": "../../../source/functions/ParsedOptimizationFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = 'r0 * sin((atan(y/x)))'
    param_symbol_names<<<{"description": "Names of parameters in 'expression' being optimized."}>>> = 'r0 r1'
    param_vector_name<<<{"description": "Reporter or VectorPostprocessor vector containing parameter values."}>>> = 'params/radii'
  []
  [h]
    type = ParsedOptimizationFunction<<<{"description": "Function used for optimization that uses a parsed expression with parameter dependence.", "href": "../../../source/functions/ParsedOptimizationFunction.html"}>>>
    # r0+${inner_radius} is the true current inner radius
    expression<<<{"description": "The user defined function."}>>> = '10 /(pi * (r0+${inner_radius})^3)'
    param_symbol_names<<<{"description": "Names of parameters in 'expression' being optimized."}>>> = 'r0 r1'
    param_vector_name<<<{"description": "Reporter or VectorPostprocessor vector containing parameter values."}>>> = 'params/radii'
  []
  [eq_grad_r0]
    type = ParsedOptimizationFunction<<<{"description": "Function used for optimization that uses a parsed expression with parameter dependence.", "href": "../../../source/functions/ParsedOptimizationFunction.html"}>>>
    expression<<<{"description": "The user defined function."}>>> = '-2 * pi * (r0 + ${inner_radius})'
    param_symbol_names<<<{"description": "Names of parameters in 'expression' being optimized."}>>> = 'r0 r1'
    param_vector_name<<<{"description": "Reporter or VectorPostprocessor vector containing parameter values."}>>> = 'params/radii'
  []
  [eq_grad_r1]
    type = ParsedOptimizationFunction<<<{"description": "Function used for optimization that uses a parsed expression with parameter dependence.", "href": "../../../source/functions/ParsedOptimizationFunction.html"}>>>
    # r1+${outer_radius} is the true current outer radius
    expression<<<{"description": "The user defined function."}>>> = '2 * pi * (r1+${outer_radius})'
    param_symbol_names<<<{"description": "Names of parameters in 'expression' being optimized."}>>> = 'r0 r1'
    param_vector_name<<<{"description": "Reporter or VectorPostprocessor vector containing parameter values."}>>> = 'params/radii'
  []
[]

[BCs<<<{"href": "../../../syntax/BCs/index.html"}>>>]
  [diffuse_r1_x]
    type = ADFunctionDirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is calculated by a function.", "href": "../../../source/bcs/ADFunctionDirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'outer'
    function<<<{"description": "The function describing the Dirichlet condition"}>>> = r1_x
    preset<<<{"description": "Whether or not to preset the BC (apply the value before the solve begins)."}>>> = false
  []
  [diffuse_r1_y]
    type = ADFunctionDirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is calculated by a function.", "href": "../../../source/bcs/ADFunctionDirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'outer'
    function<<<{"description": "The function describing the Dirichlet condition"}>>> = r1_y
    preset<<<{"description": "Whether or not to preset the BC (apply the value before the solve begins)."}>>> = false
  []
  [diffuse_r0_x]
    type = ADFunctionDirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is calculated by a function.", "href": "../../../source/bcs/ADFunctionDirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'inner'
    function<<<{"description": "The function describing the Dirichlet condition"}>>> = r0_x
    preset<<<{"description": "Whether or not to preset the BC (apply the value before the solve begins)."}>>> = false
  []
  [diffuse_r0_y]
    type = ADFunctionDirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is calculated by a function.", "href": "../../../source/bcs/ADFunctionDirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = 'inner'
    function<<<{"description": "The function describing the Dirichlet condition"}>>> = r0_y
    preset<<<{"description": "Whether or not to preset the BC (apply the value before the solve begins)."}>>> = false
  []
  # run physics on deformed mesh
  [convection]
    type = ADMatNeumannBC<<<{"description": "Imposes the integrated boundary condition $\\frac{C \\partial u}{\\partial n}=M*h$, where $h$ is a constant, $M$ is a material property, and $C$ is a coefficient defined by the kernel for $u$.", "href": "../../../source/bcs/MatNeumannBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = inner
    boundary_material<<<{"description": "Material property multiplying the constant that will be enforced by the BC"}>>> = convection
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation.  Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = true
    value<<<{"description": "For a Laplacian problem, the value of the gradient dotted with the normals on the boundary."}>>> = 1
  []
[]

[Kernels<<<{"href": "../../../syntax/Kernels/index.html"}>>>]
  [disp_x]
    type = MatDiffusion<<<{"description": "Diffusion equation Kernel that takes an isotropic Diffusivity from a material property", "href": "../../../source/kernels/MatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_x
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = diff_coef
  []
  [disp_y]
    type = MatDiffusion<<<{"description": "Diffusion equation Kernel that takes an isotropic Diffusivity from a material property", "href": "../../../source/kernels/MatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = disp_y
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = false
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = diff_coef
  []
  # run physics of interest on deformed mesh
  [Diffusion]
    type = FunctionDiffusion<<<{"description": "Diffusion with a function coefficient.", "href": "../../../source/kernels/FunctionDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = T
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = true
  []
  [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"}>>> = T
    value<<<{"description": "Coefficient to multiply by the body force term"}>>> = 1
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = true
  []
[]

For the equality constraints, a gradient of the constraint with respect to the parameters is needed for TAO. The gradient is analytically computed for the current radii.

Using postprocessors, the objective and constraint functions can be calculated. The objective in this case is to minimize the maximum temperature. Also the current volume is calculated and how much that volume violates the constraint. These values are now able to be transferred to the main application for the optimization process.

Listing 2: Forward App Objective Calculations

[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [current_volume]
    type = VolumePostprocessor<<<{"description": "Computes the volume of a specified block", "href": "../../../source/postprocessors/VolumePostprocessor.html"}>>>
    use_displaced_mesh<<<{"description": "Whether or not this object should use the displaced mesh for computation.  Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used."}>>> = true
  []
  # objective function
  [objective]
    type = NodalExtremeValue<<<{"description": "Finds either the min or max elemental value of a variable over the domain.", "href": "../../../source/postprocessors/NodalExtremeValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on"}>>> = T
  []
  [eq_constraint]
    type = ParsedPostprocessor<<<{"description": "Computes a parsed expression with post-processors", "href": "../../../source/postprocessors/ParsedPostprocessor.html"}>>>
    pp_names<<<{"description": "Post-processors arguments"}>>> = current_volume
    expression<<<{"description": "function expression"}>>> = 'current_volume - ${volume_constraint}'
  []
  [func_r0]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = eq_grad_r0
  []
  [func_r1]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = eq_grad_r1
  []
[]

Main-App Optimization

Optimization problems are solved using the MultiApps system. The main application contains the optimization executioner and the sub-application solve the forward problem.

For this optimization example, the GeneralOptimization reporter is used. The "equality_names" option lists equality constraints. There is also a "inequality_names", not shown below, for listing the inequality constraints. The GeneralOptimization also expects the controllable parameters "objective_name" and "num_values_name", which are the names of the reporters that will hold the objective value and the number of parameters in the forward problem.

Listing 3: Main application OptimizationReporter block with equality constraints

[OptimizationReporter<<<{"href": "../../../syntax/OptimizationReporter/index.html"}>>>]
  type = GeneralOptimization
  parameter_names = 'radii'
  num_values_name = num_radii
  equality_names = 'volume_constraint'
  initial_condition = '0 0'
  lower_bounds = '-10'
  upper_bounds = '10'
  objective_name = max_temp
[]

To enable constrained optimization the "tao_solver" needs to be set to taoalmm, which is an Augmented Lagrangian method. This is the only TAO solver in MOOSE that has constraint support. The taoalmm subsolver is set with "petsc_options_iname" and "petsc_options_value". This example computes the objective function's gradient using finite differencing in TAO with the petsc_options_iname -tao_fd. Finite differencing only uses the forward problem's objective value and does not require an adjoint solve. -tao_almm_subsolver_tao_type and bqnktr set the subsolver to the Bounded Quasi-Newton-Krylov Trust Region Method.

Listing 4: Main application Executioner block

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Optimize
  tao_solver = taoalmm
  petsc_options_iname = '-tao_almm_subsolver_tao_type -tao_gatol -tao_catol  -tao_almm_type -tao_almm_mu_init -tao_fd_gradient -tao_fd_delta'
  petsc_options_value = 'bqnktr 1e-3 1e-3  phr 1e9 true 1e-8 '
  verbose = true
[]

Adjoint Problem Sub-App

Since we are using finite-differencing to determine the gradients of our objective function, a Adjoint Sub-App is not needed.

Optimization Results

As shown in Figure 1, the initial annulus design has a temperature profile that reaches a maximum temperature of around 470. The initial inner and outer radii for this design were 6 and 10, respectively.

Figure 1: Initial annulus temperature profile

After optimizing the design to decrease the maximum temperature, the new design has a maximum temperature of around 160. Moreover, the volume constraint was satisfied to the given tolerance. The optimized inner and outer radii are now 1.199 and 8.075, respectively.

Figure 2 shows the optimized designs and the new temperature profile.

Figure 2: Optimized Annulus temperature profile

(modules/optimization/test/tests/executioners/constrained/shape_optimization/forward.i)

# This test is documented as an example for ConstrainedShapeOptimization. This
# test should not be changed without updating the documentation.

inner_radius = 6
outer_radius = 10
volume_constraint = 200

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [mesh]
    type = ConcentricCircleMeshGenerator
    has_outer_square = no
    num_sectors = 16
    radii = '${inner_radius} ${outer_radius}'
    rings = '16 16'
    preserve_volumes = false
  []
  [inner_radius]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    new_boundary = inner
    primary_block = 2
    paired_block = 1
  []

  [delete]
    type = BlockDeletionGenerator
    input = inner_radius
    block = 1
  []
  [gather_all]
    type = BoundingBoxNodeSetGenerator
    input = delete
    bottom_left = '-100 -100 -100'
    top_right = '100 100 100'
    new_boundary = total
  []
  [combine]
    type = SideSetsFromBoundingBoxGenerator
    input = gather_all
    bottom_left = '-100 -100 -100'
    top_right = '100 100 100'
    boundaries_old = 'inner outer'
    boundary_new = moving
  []

[]
[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [T]
  []
[]

[AuxVariables]
  [dist_between]
    [AuxKernel]
      type = NearestNodeDistanceAux
      variable = dist_between
      paired_boundary = moving
      boundary = total
      block = 2
      use_displaced_mesh = false
      execute_on = "INITIAL"
    []
  []
[]

[Kernels]
  [disp_x]
    type = MatDiffusion
    variable = disp_x
    use_displaced_mesh = false
    diffusivity = diff_coef
  []
  [disp_y]
    type = MatDiffusion
    variable = disp_y
    use_displaced_mesh = false
    diffusivity = diff_coef
  []
  # run physics of interest on deformed mesh
  [Diffusion]
    type = FunctionDiffusion
    variable = T
    use_displaced_mesh = true
  []
  [Source]
    type = BodyForce
    variable = T
    value = 1
    use_displaced_mesh = true
  []
[]

[Materials]
  # perserve elements near the boundary
  [diff_coef]
    type = ParsedMaterial
    coupled_variables = 'dist_between'
    expression = '1/(dist_between+0.5)'
    property_name = 'diff_coef'
  []
  [h]
    type = ADGenericFunctionMaterial
    prop_names = h
    prop_values = h
  []
  # convection type boundary condition
  [convection_bc]
    type = ADParsedMaterial
    coupled_variables = "T"
    expression = "h*(100-T)"
    material_property_names = "h"
    property_name = convection
  []
[]

[Functions]
  [r1_x]
    type = ParsedOptimizationFunction
    expression = 'r1 * cos((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [r1_y]
    type = ParsedOptimizationFunction
    expression = 'r1 * sin((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [r0_x]
    type = ParsedOptimizationFunction
    expression = 'r0 * cos((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [r0_y]
    type = ParsedOptimizationFunction
    expression = 'r0 * sin((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [h]
    type = ParsedOptimizationFunction
    # r0+${inner_radius} is the true current inner radius
    expression = '10 /(pi * (r0+${inner_radius})^3)'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [eq_grad_r0]
    type = ParsedOptimizationFunction
    expression = '-2 * pi * (r0 + ${inner_radius})'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [eq_grad_r1]
    type = ParsedOptimizationFunction
    # r1+${outer_radius} is the true current outer radius
    expression = '2 * pi * (r1+${outer_radius})'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
[]

[BCs]
  [diffuse_r1_x]
    type = ADFunctionDirichletBC
    variable = disp_x
    boundary = 'outer'
    function = r1_x
    preset = false
  []
  [diffuse_r1_y]
    type = ADFunctionDirichletBC
    variable = disp_y
    boundary = 'outer'
    function = r1_y
    preset = false
  []
  [diffuse_r0_x]
    type = ADFunctionDirichletBC
    variable = disp_x
    boundary = 'inner'
    function = r0_x
    preset = false
  []
  [diffuse_r0_y]
    type = ADFunctionDirichletBC
    variable = disp_y
    boundary = 'inner'
    function = r0_y
    preset = false
  []
  # run physics on deformed mesh
  [convection]
    type = ADMatNeumannBC
    variable = T
    boundary = inner
    boundary_material = convection
    use_displaced_mesh = true
    value = 1
  []
[]

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

[Reporters]
  [params]
    type = ConstantReporter
    real_vector_names = 'radii'
    real_vector_values = '0 0'
    dof_id_type_vector_names = 'num_params'
    dof_id_type_vector_values = '2'
  []
[]

[Postprocessors]
  [current_volume]
    type = VolumePostprocessor
    use_displaced_mesh = true
  []
  # objective function
  [objective]
    type = NodalExtremeValue
    variable = T
  []
  [eq_constraint]
    type = ParsedPostprocessor
    pp_names = current_volume
    expression = 'current_volume - ${volume_constraint}'
  []
  [func_r0]
    type = FunctionValuePostprocessor
    function = eq_grad_r0
  []
  [func_r1]
    type = FunctionValuePostprocessor
    function = eq_grad_r1
  []
[]

[VectorPostprocessors]
  # convert "Real" postprocessors to vectors
  [vol_constraint]
    type = VectorOfPostprocessors
    postprocessors = 'eq_constraint'
    force_postaux = true
  []
  [eq_grad]
    type = VectorOfPostprocessors
    postprocessors = 'func_r0 func_r1'
    force_postaux = true
  []
[]

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  line_search = none
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-8
[]

[Outputs]
  console = false
[]

(modules/optimization/test/tests/executioners/constrained/shape_optimization/forward.i)

# This test is documented as an example for ConstrainedShapeOptimization. This
# test should not be changed without updating the documentation.

inner_radius = 6
outer_radius = 10
volume_constraint = 200

[GlobalParams]
  displacements = 'disp_x disp_y'
[]

[Mesh]
  [mesh]
    type = ConcentricCircleMeshGenerator
    has_outer_square = no
    num_sectors = 16
    radii = '${inner_radius} ${outer_radius}'
    rings = '16 16'
    preserve_volumes = false
  []
  [inner_radius]
    type = SideSetsBetweenSubdomainsGenerator
    input = mesh
    new_boundary = inner
    primary_block = 2
    paired_block = 1
  []

  [delete]
    type = BlockDeletionGenerator
    input = inner_radius
    block = 1
  []
  [gather_all]
    type = BoundingBoxNodeSetGenerator
    input = delete
    bottom_left = '-100 -100 -100'
    top_right = '100 100 100'
    new_boundary = total
  []
  [combine]
    type = SideSetsFromBoundingBoxGenerator
    input = gather_all
    bottom_left = '-100 -100 -100'
    top_right = '100 100 100'
    boundaries_old = 'inner outer'
    boundary_new = moving
  []

[]
[Variables]
  [disp_x]
  []
  [disp_y]
  []
  [T]
  []
[]

[AuxVariables]
  [dist_between]
    [AuxKernel]
      type = NearestNodeDistanceAux
      variable = dist_between
      paired_boundary = moving
      boundary = total
      block = 2
      use_displaced_mesh = false
      execute_on = "INITIAL"
    []
  []
[]

[Kernels]
  [disp_x]
    type = MatDiffusion
    variable = disp_x
    use_displaced_mesh = false
    diffusivity = diff_coef
  []
  [disp_y]
    type = MatDiffusion
    variable = disp_y
    use_displaced_mesh = false
    diffusivity = diff_coef
  []
  # run physics of interest on deformed mesh
  [Diffusion]
    type = FunctionDiffusion
    variable = T
    use_displaced_mesh = true
  []
  [Source]
    type = BodyForce
    variable = T
    value = 1
    use_displaced_mesh = true
  []
[]

[Materials]
  # perserve elements near the boundary
  [diff_coef]
    type = ParsedMaterial
    coupled_variables = 'dist_between'
    expression = '1/(dist_between+0.5)'
    property_name = 'diff_coef'
  []
  [h]
    type = ADGenericFunctionMaterial
    prop_names = h
    prop_values = h
  []
  # convection type boundary condition
  [convection_bc]
    type = ADParsedMaterial
    coupled_variables = "T"
    expression = "h*(100-T)"
    material_property_names = "h"
    property_name = convection
  []
[]

[Functions]
  [r1_x]
    type = ParsedOptimizationFunction
    expression = 'r1 * cos((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [r1_y]
    type = ParsedOptimizationFunction
    expression = 'r1 * sin((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [r0_x]
    type = ParsedOptimizationFunction
    expression = 'r0 * cos((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [r0_y]
    type = ParsedOptimizationFunction
    expression = 'r0 * sin((atan(y/x)))'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [h]
    type = ParsedOptimizationFunction
    # r0+${inner_radius} is the true current inner radius
    expression = '10 /(pi * (r0+${inner_radius})^3)'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [eq_grad_r0]
    type = ParsedOptimizationFunction
    expression = '-2 * pi * (r0 + ${inner_radius})'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
  [eq_grad_r1]
    type = ParsedOptimizationFunction
    # r1+${outer_radius} is the true current outer radius
    expression = '2 * pi * (r1+${outer_radius})'
    param_symbol_names = 'r0 r1'
    param_vector_name = 'params/radii'
  []
[]

[BCs]
  [diffuse_r1_x]
    type = ADFunctionDirichletBC
    variable = disp_x
    boundary = 'outer'
    function = r1_x
    preset = false
  []
  [diffuse_r1_y]
    type = ADFunctionDirichletBC
    variable = disp_y
    boundary = 'outer'
    function = r1_y
    preset = false
  []
  [diffuse_r0_x]
    type = ADFunctionDirichletBC
    variable = disp_x
    boundary = 'inner'
    function = r0_x
    preset = false
  []
  [diffuse_r0_y]
    type = ADFunctionDirichletBC
    variable = disp_y
    boundary = 'inner'
    function = r0_y
    preset = false
  []
  # run physics on deformed mesh
  [convection]
    type = ADMatNeumannBC
    variable = T
    boundary = inner
    boundary_material = convection
    use_displaced_mesh = true
    value = 1
  []
[]

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

[Reporters]
  [params]
    type = ConstantReporter
    real_vector_names = 'radii'
    real_vector_values = '0 0'
    dof_id_type_vector_names = 'num_params'
    dof_id_type_vector_values = '2'
  []
[]

[Postprocessors]
  [current_volume]
    type = VolumePostprocessor
    use_displaced_mesh = true
  []
  # objective function
  [objective]
    type = NodalExtremeValue
    variable = T
  []
  [eq_constraint]
    type = ParsedPostprocessor
    pp_names = current_volume
    expression = 'current_volume - ${volume_constraint}'
  []
  [func_r0]
    type = FunctionValuePostprocessor
    function = eq_grad_r0
  []
  [func_r1]
    type = FunctionValuePostprocessor
    function = eq_grad_r1
  []
[]

[VectorPostprocessors]
  # convert "Real" postprocessors to vectors
  [vol_constraint]
    type = VectorOfPostprocessors
    postprocessors = 'eq_constraint'
    force_postaux = true
  []
  [eq_grad]
    type = VectorOfPostprocessors
    postprocessors = 'func_r0 func_r1'
    force_postaux = true
  []
[]

[Executioner]
  type = Steady
  solve_type = 'NEWTON'
  petsc_options_iname = '-pc_type -pc_factor_shift_type'
  petsc_options_value = 'lu       NONZERO'
  line_search = none
  nl_rel_tol = 1e-8
  nl_abs_tol = 1e-8
[]

[Outputs]
  console = false
[]

(modules/optimization/test/tests/executioners/constrained/shape_optimization/main.i)

[Optimization]
[]

[OptimizationReporter]
  type = GeneralOptimization
  parameter_names = 'radii'
  num_values_name = num_radii
  equality_names = 'volume_constraint'
  initial_condition = '0 0'
  lower_bounds = '-10'
  upper_bounds = '10'
  objective_name = max_temp
[]

[Executioner]
  type = Optimize
  tao_solver = taoalmm
  petsc_options_iname = '-tao_almm_subsolver_tao_type -tao_gatol -tao_catol  -tao_almm_type -tao_almm_mu_init -tao_fd_gradient -tao_fd_delta'
  petsc_options_value = 'bqnktr 1e-3 1e-3  phr 1e9 true 1e-8 '
  verbose = true
[]

[MultiApps]
  [forward]
    type = FullSolveMultiApp
    input_files = forward.i
    # Run on initial so the forward problem can determine number of parameters
    execute_on = 'INITIAL FORWARD'
  []
[]

[Transfers]
  [to_forward]
    type = MultiAppReporterTransfer
    to_multi_app = forward
    from_reporters = 'OptimizationReporter/radii'
    to_reporters = 'params/radii'
    execute_on = FORWARD
  []
  [from_forward]
    type = MultiAppReporterTransfer
    from_multi_app = forward
    from_reporters = 'params/num_params
                      objective/value
                      vol_constraint/vol_constraint
                      eq_grad/eq_grad'
    to_reporters = 'OptimizationReporter/num_radii
                    OptimizationReporter/max_temp
                    OptimizationReporter/volume_constraint
                    OptimizationReporter/grad_volume_constraint'
  []
[]

[Outputs]
  csv = true
  print_linear_residuals = false
  print_nonlinear_residuals = false
[]

(modules/optimization/test/tests/executioners/constrained/shape_optimization/main.i)

[Optimization]
[]

[OptimizationReporter]
  type = GeneralOptimization
  parameter_names = 'radii'
  num_values_name = num_radii
  equality_names = 'volume_constraint'
  initial_condition = '0 0'
  lower_bounds = '-10'
  upper_bounds = '10'
  objective_name = max_temp
[]

[Executioner]
  type = Optimize
  tao_solver = taoalmm
  petsc_options_iname = '-tao_almm_subsolver_tao_type -tao_gatol -tao_catol  -tao_almm_type -tao_almm_mu_init -tao_fd_gradient -tao_fd_delta'
  petsc_options_value = 'bqnktr 1e-3 1e-3  phr 1e9 true 1e-8 '
  verbose = true
[]

[MultiApps]
  [forward]
    type = FullSolveMultiApp
    input_files = forward.i
    # Run on initial so the forward problem can determine number of parameters
    execute_on = 'INITIAL FORWARD'
  []
[]

[Transfers]
  [to_forward]
    type = MultiAppReporterTransfer
    to_multi_app = forward
    from_reporters = 'OptimizationReporter/radii'
    to_reporters = 'params/radii'
    execute_on = FORWARD
  []
  [from_forward]
    type = MultiAppReporterTransfer
    from_multi_app = forward
    from_reporters = 'params/num_params
                      objective/value
                      vol_constraint/vol_constraint
                      eq_grad/eq_grad'
    to_reporters = 'OptimizationReporter/num_radii
                    OptimizationReporter/max_temp
                    OptimizationReporter/volume_constraint
                    OptimizationReporter/grad_volume_constraint'
  []
[]

[Outputs]
  csv = true
  print_linear_residuals = false
  print_nonlinear_residuals = false
[]

Problem Description
Forward Problem Sub - App
Main - App Optimization
Adjoint Problem Sub - App
Optimization Results