Force Inversion Example: Neumann Boundary Condition

Background

The MOOSE optimization module provides a flexible framework for solving inverse optimization problems in MOOSE. This page is part of a set of examples for different types of inverse optimization problems.

Example: Neumann Boundary Condition

Main-App Optimization Executioner

In this example, we will be optimizing the multi-objective function given by:

(1)

where represents the displacement vector, denotes the parameters, and are the observed and target -components of the displacement at the -th point, and and are the observed and target -components of the displacement at the -th point.

The derivative of the objective function with respect to the parameters can be expressed as:

(2)

In this example we will be optimizing the two forcing parameters on the right side of a cantilever beam. The forcing parameters control the force that is applied to the beam in each direction.

The main app, shown in Listing 1, only receives the final objective and gradient.

Listing 1: MainApp Input File

[Optimization<<<{"href": "../../../syntax/Optimization/index.html"}>>>]
[]

[Mesh<<<{"href": "../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  nx = 10
  ny = 2
  xmin = 0.0
  xmax = 5.0
  ymin = 0.0
  ymax = 1.0
[]

[OptimizationReporter<<<{"href": "../../../syntax/OptimizationReporter/index.html"}>>>]
  type = GeneralOptimization
  objective_name = objective_value
  parameter_names = 'force_right'
  num_values = '2'
  initial_condition = '100'
  outputs = 'csv'
[]

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Optimize
  tao_solver = taobqnktr
  petsc_options_iname = '-tao_gttol -tao_max_it'
  petsc_options_value = '1e-7 500 '
  verbose = true
[]

[MultiApps<<<{"href": "../../../syntax/MultiApps/index.html"}>>>]
  [forward]
    type = FullSolveMultiApp<<<{"description": "Performs a complete simulation during each execution.", "href": "../../../source/multiapps/FullSolveMultiApp.html"}>>>
    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."}>>> = forward_and_adjoint.i
    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."}>>> = "FORWARD"
    clone_parent_mesh<<<{"description": "True to clone parent app mesh and use it for this MultiApp."}>>> = true
  []
[]

[Transfers<<<{"href": "../../../syntax/Transfers/index.html"}>>>]
  [toForward_measument]
    type = MultiAppReporterTransfer<<<{"description": "Transfers reporter data between two applications.", "href": "../../../source/transfers/MultiAppReporterTransfer.html"}>>>
    to_multi_app<<<{"description": "The name of the MultiApp to transfer the data to"}>>> = forward
    from_reporters<<<{"description": "List of the reporter names (object_name/value_name) to transfer the value from."}>>> = 'OptimizationReporter/force_right'
    to_reporters<<<{"description": "List of the reporter names (object_name/value_name) to transfer the value to."}>>> = 'params/right_values'
  []
  [fromForward]
    type = MultiAppReporterTransfer<<<{"description": "Transfers reporter data between two applications.", "href": "../../../source/transfers/MultiAppReporterTransfer.html"}>>>
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = forward
    from_reporters<<<{"description": "List of the reporter names (object_name/value_name) to transfer the value from."}>>> = 'combined/gradient obj/obj_val'
    to_reporters<<<{"description": "List of the reporter names (object_name/value_name) to transfer the value to."}>>> = 'OptimizationReporter/grad_force_right OptimizationReporter/objective_value'
  []
[]

[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'
[]
(modules/combined/test/tests/optimization/invOpt_mechanics/main_auto_adjoint.i)

Forward and Adjoint Problem Sub-App

The definition of the objective function in Eq. (1) allows for the gradient to be the addition of two gradients: one for the -component and one for the -component. The forward sub app is responsible for correctly creating both sub-gradients and combining them, as well as the objective functions, to be used by the main app. This example will only go into detail on how to calculate and combine these objectives and gradients. For more details on setting up the entire optimization process, look at Example 1: Point Loads.

For functions, we split the parameter vector, which has two entries, into two separate functions: one that controls the force for the -component and another that controls the force for the -component.

Listing 2: Functions

[Functions<<<{"href": "../../../syntax/Functions/index.html"}>>>]
  [right_fy_func]
    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."}>>> = 'val_y'
    param_symbol_names<<<{"description": "Names of parameters in 'expression' being optimized."}>>> = 'val_x val_y'
    param_vector_name<<<{"description": "Reporter or VectorPostprocessor vector containing parameter values."}>>> = 'params/right_values'
  []
  [right_fx_func]
    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."}>>> = 'val_x'
    param_symbol_names<<<{"description": "Names of parameters in 'expression' being optimized."}>>> = 'val_x val_y'
    param_vector_name<<<{"description": "Reporter or VectorPostprocessor vector containing parameter values."}>>> = 'params/right_values'
  []
[]
(modules/combined/test/tests/optimization/invOpt_mechanics/forward_and_adjoint.i)

For the boundary conditions, we apply the forces as Neumann boundary conditions.

Listing 3: Boundary Conditions

[BCs<<<{"href": "../../../syntax/BCs/index.html"}>>>]
  [left_ux]
    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"}>>> = disp_x
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 0
  []
  [left_uy]
    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"}>>> = disp_y
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = 0
  []
  [right_fy]
    type = FunctionNeumannBC<<<{"description": "Imposes the integrated boundary condition $\\frac{\\partial u}{\\partial n}=h(t,\\vec{x})$, where $h$ is a (possibly) time and space-dependent MOOSE Function.", "href": "../../../source/bcs/FunctionNeumannBC.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"}>>> = right
    function<<<{"description": "The function."}>>> = right_fy_func
  []
  [right_fx]
    type = FunctionNeumannBC<<<{"description": "Imposes the integrated boundary condition $\\frac{\\partial u}{\\partial n}=h(t,\\vec{x})$, where $h$ is a (possibly) time and space-dependent MOOSE Function.", "href": "../../../source/bcs/FunctionNeumannBC.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"}>>> = right
    function<<<{"description": "The function."}>>> = right_fx_func
  []
[]
(modules/combined/test/tests/optimization/invOpt_mechanics/forward_and_adjoint.i)

For the reporters, we create two OptimizationData objects that hold the observed measurement data for and individually. This Reporter will calculate the individual objective values. The locations or the number of points for and don't need to match. Also, we use a ParsedVectorReporter and a ParsedScalarReporter to combine the individual gradients and objective function, respectively. These are the final reporters that will be transferred to the main app.

Listing 4: Reporters

[Reporters<<<{"href": "../../../syntax/Reporters/index.html"}>>>]
  [measure_data_x]
    type = OptimizationData<<<{"description": "Reporter to hold measurement and simulation data for optimization problems", "href": "../../../source/reporters/OptimizationData.html"}>>>
    objective_name<<<{"description": "Name of reporter value defining the objective."}>>> = objective_value
    variable<<<{"description": "Vector of variable names to sample at measurement points."}>>> = disp_x
    measurement_points<<<{"description": "Point locations corresponding to each measurement value"}>>> = '5.0 1.0 0.0'
    measurement_values<<<{"description": "Measurement values collected from locations given by measurement_points"}>>> = '-13.00873469088'
  []
  [measure_data_y]
    type = OptimizationData<<<{"description": "Reporter to hold measurement and simulation data for optimization problems", "href": "../../../source/reporters/OptimizationData.html"}>>>
    objective_name<<<{"description": "Name of reporter value defining the objective."}>>> = objective_value
    variable<<<{"description": "Vector of variable names to sample at measurement points."}>>> = disp_y
    measurement_points<<<{"description": "Point locations corresponding to each measurement value"}>>> = '5.0 1.0 0.0'
    measurement_values<<<{"description": "Measurement values collected from locations given by measurement_points"}>>> = '85.008027719915'
  []
  [params]
    type = ConstantReporter<<<{"description": "Reporter with constant values to be accessed by other objects, can be modified using transfers.", "href": "../../../source/reporters/ConstantReporter.html"}>>>
    real_vector_names<<<{"description": "Names for each vector of reals value."}>>> = 'right_values'
    real_vector_values<<<{"description": "Values for vectors of reals."}>>> = '-1300 2100 ' # True Values
  []
  [combined]
    type = ParsedVectorReporter<<<{"description": "Apply parsed functions to vector entries held in reporters.", "href": "../../../source/reporters/ParsedVectorReporter.html"}>>>
    name<<<{"description": "Name of output reporter."}>>> = gradient
    reporter_names<<<{"description": "Reporter names "}>>> = 'adjoint_pt_x/inner_product adjoint_pt_y/inner_product'
    reporter_symbols<<<{"description": "Expression symbol for each reporter"}>>> = 'a b'
    expression<<<{"description": "function expression"}>>> = 'a+b'
    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."}>>> = ADJOINT_TIMESTEP_END
    # Just to confirm this happens after the gradient calcutions
    execution_order_group<<<{"description": "Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group."}>>> = 1
  []
  [obj]
    type = ParsedScalarReporter<<<{"description": "Applies parsed functions to scalar entries held in reporters.", "href": "../../../source/reporters/ParsedScalarReporter.html"}>>>
    name<<<{"description": "Name of output reporter."}>>> = obj_val
    reporter_names<<<{"description": "The input reporter names"}>>> = 'measure_data_x/objective_value
    measure_data_y/objective_value'
    reporter_symbols<<<{"description": "Expression symbol for each reporter"}>>> = 'a b'
    expression<<<{"description": "function expression"}>>> = 'a+b'
    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."}>>> = ADJOINT_TIMESTEP_END
    # Just to confirm this happens after the gradient calcutions
    execution_order_group<<<{"description": "Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group."}>>> = 1
  []
[]
(modules/combined/test/tests/optimization/invOpt_mechanics/forward_and_adjoint.i)

The Dirac kernels show how we apply the misfit to the adjoint variables, with the adjoint getting the misfit from the data and the getting the misfit from the data.

Listing 5: Dirac Kernels

[DiracKernels<<<{"href": "../../../syntax/DiracKernels/index.html"}>>>]
  [pt_x]
    type = ReporterPointSource<<<{"description": "Apply a point load defined by Reporter.", "href": "../../../source/dirackernels/ReporterPointSource.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = adjoint_x
    x_coord_name<<<{"description": "reporter x-coordinate name.  This uses the reporter syntax <reporter>/<name>."}>>> = measure_data_x/measurement_xcoord
    y_coord_name<<<{"description": "reporter y-coordinate name.  This uses the reporter syntax <reporter>/<name>."}>>> = measure_data_x/measurement_ycoord
    z_coord_name<<<{"description": "reporter z-coordinate name.  This uses the reporter syntax <reporter>/<name>."}>>> = measure_data_x/measurement_zcoord
    value_name<<<{"description": "reporter value name.  This uses the reporter syntax <reporter>/<name>."}>>> = measure_data_x/misfit_values
  []
  [pt_y]
    type = ReporterPointSource<<<{"description": "Apply a point load defined by Reporter.", "href": "../../../source/dirackernels/ReporterPointSource.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = adjoint_y
    x_coord_name<<<{"description": "reporter x-coordinate name.  This uses the reporter syntax <reporter>/<name>."}>>> = measure_data_y/measurement_xcoord
    y_coord_name<<<{"description": "reporter y-coordinate name.  This uses the reporter syntax <reporter>/<name>."}>>> = measure_data_y/measurement_ycoord
    z_coord_name<<<{"description": "reporter z-coordinate name.  This uses the reporter syntax <reporter>/<name>."}>>> = measure_data_y/measurement_zcoord
    value_name<<<{"description": "reporter value name.  This uses the reporter syntax <reporter>/<name>."}>>> = measure_data_y/misfit_values
  []
[]
(modules/combined/test/tests/optimization/invOpt_mechanics/forward_and_adjoint.i)

Finally, in vector postprocessors, we show that we calculate the gradients for both components using SideOptimizationNeumannFunctionInnerProduct with the correct adjoint variable and function. These gradients are then combined in the reporters block as previously discussed.

Listing 6: Vector Postprocessors

[VectorPostprocessors<<<{"href": "../../../syntax/VectorPostprocessors/index.html"}>>>]
  [point_sample]
    type = PointValueSampler<<<{"description": "Sample a variable at specific points.", "href": "../../../source/vectorpostprocessors/PointValueSampler.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = 'disp_x disp_y'
    points<<<{"description": "The points where you want to evaluate the variables"}>>> = '5.0 1.0 0'
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
  []

  [adjoint_pt_x]
    type = SideOptimizationNeumannFunctionInnerProduct<<<{"description": "Computes the inner product of variable with parameterized Neumann function for optimization gradient computation.", "href": "../../../source/vectorpostprocessors/SideOptimizationNeumannFunctionInnerProduct.html"}>>>
    variable<<<{"description": "Variable used for inner product calculation, usually the adjoint variable."}>>> = adjoint_x
    function<<<{"description": "Optimization function."}>>> = right_fx_func
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    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."}>>> = ADJOINT_TIMESTEP_END
  []
  [adjoint_pt_y]
    type = SideOptimizationNeumannFunctionInnerProduct<<<{"description": "Computes the inner product of variable with parameterized Neumann function for optimization gradient computation.", "href": "../../../source/vectorpostprocessors/SideOptimizationNeumannFunctionInnerProduct.html"}>>>
    variable<<<{"description": "Variable used for inner product calculation, usually the adjoint variable."}>>> = adjoint_y
    function<<<{"description": "Optimization function."}>>> = right_fy_func
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    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."}>>> = ADJOINT_TIMESTEP_END
  []
[]
(modules/combined/test/tests/optimization/invOpt_mechanics/forward_and_adjoint.i)