Material Inversion Example: Nonlinear Diffusion Reaction

The following illustrates the capability of the MOOSE Optimization module by applying the module to a nonlinear, material-inversion optimization problem:

var element = document.getElementById("moose-equation-e72cef1b-7cab-4ef8-8353-a2ff01382a7f");katex.render("\\begin{split} &\\min_p f(u, \\sigma) = \\frac{1}{2}\\sum^N_{i=1} \\left(u_i - \\tilde{u}_i\\right)^2; \\\\ &\\text{subject to} \\left\\{\\begin{array}{lr} \\dot{u} - \\vec{\\nabla}\\cdot(1+u)\\vec{\\nabla}u + \\sigma u = 1 & \\text{in} \\, t=[0,1],\\Omega=(0,1)^2\\\\ u(t=0) = 0 & \\\\ u(x=0) = u(y=0) = 0 & \\\\ \\left.\\frac{\\partial u}{\\partial x}\\right|_{x=1} = \\left.\\frac{\\partial u}{\\partial y}\\right|_{y=1} = 0 \\\\ \\sigma \\geq 0 \\end{array}\\right., \\end{split}", 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 $var element = document.getElementById("moose-equation-c73983b7-571c-4b7b-8b14-a1cef4f08557");katex.render("N", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ is the number of measurement locations and $var element = document.getElementById("moose-equation-b825a5be-12ad-457c-bb72-c3931cff6cfd");katex.render("u_i", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and $var element = document.getElementById("moose-equation-09bd4462-8b38-446e-b7aa-753139b77d79");katex.render("\\tilde{u}_i", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ are the simulated and measured state variables (e.g. temperature and displacement fields), respectively, at location $var element = document.getElementById("moose-equation-c936410f-9c3a-4846-bf2a-74f35330920e");katex.render("i", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ . The parameter being optimized is the spatially dependent reaction rate ( $var element = document.getElementById("moose-equation-998f74e0-d2be-41aa-a47f-50339a8da075");katex.render("\\sigma", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ ). The domain is meshed using a $var element = document.getElementById("moose-equation-30e045f1-059a-45aa-bfa1-d64aa36d87ea");katex.render("16\\times16", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ grid of quadrilateral elements, and time is discretized using implict Euler over ten uniform time steps. The measurement data is generated by evaluating the PDE with $var element = document.getElementById("moose-equation-2b33ffe6-57c1-48d0-9862-fb49e664f934");katex.render("\\sigma = e^{xy} - 1", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ and sampling the resulting solution at 22 locations—shown in top left plot of Figure 1—at every time step, resulting in $var element = document.getElementById("moose-equation-033fa3b9-d5c0-49a4-8127-031b3373e636");katex.render("N=220", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ .

The reaction rate is parameterized using the mesh shown in the top right plot of Figure 1, where parameter values set the reaction rate at the 19 nodes and linearly interpolating between them. The initial condition for the optimization sets $var element = document.getElementById("moose-equation-6247cf02-ccc9-4f14-a111-43a74be0bc6e");katex.render("\\sigma=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}"}});$ , emulating a diffusion-only system. TAO's bounded quasi-Newton line search (TAOBQNLS) was the chosen optimization algorithm. The following listings show the optimization input driving the optimization solve and the physics sub-application input used to calculate the forward and adjoint simulations.

Listing 1: Optimization input

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

[OptimizationReporter<<<{"href": "../../../syntax/OptimizationReporter/index.html"}>>>]
  type = ParameterMeshOptimization
  objective_name = objective_value
  parameter_names = 'reaction_rate'
  parameter_meshes = 'parameter_mesh_out.e'

  initial_condition = 0
  lower_bounds = 0

[]

[Reporters<<<{"href": "../../../syntax/Reporters/index.html"}>>>]
  [main]
    type = OptimizationData<<<{"description": "Reporter to hold measurement and simulation data for optimization problems", "href": "../../../source/reporters/OptimizationData.html"}>>>
    measurement_file<<<{"description": "CSV file with measurement value and coordinates (value, x, y, z)."}>>> = forward_exact_csv_sample_0011.csv
    file_xcoord<<<{"description": "x coordinate column name from measurement_file csv being read in."}>>> = measurement_xcoord
    file_ycoord<<<{"description": "y coordinate column name from csv file being read in."}>>> = measurement_ycoord
    file_zcoord<<<{"description": "z coordinate column name from csv file being read in."}>>> = measurement_zcoord
    file_time<<<{"description": "time column name from csv file being read in."}>>> = measurement_time
    file_value<<<{"description": "measurement value column name from csv file being read in."}>>> = simulation_values
  []
[]

[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
  []
[]

[Transfers<<<{"href": "../../../syntax/Transfers/index.html"}>>>]
  [to_forward]
    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."}>>> = 'main/measurement_xcoord
                      main/measurement_ycoord
                      main/measurement_zcoord
                      main/measurement_time
                      main/measurement_values
                      OptimizationReporter/reaction_rate'
    to_reporters<<<{"description": "List of the reporter names (object_name/value_name) to transfer the value to."}>>> = 'data/measurement_xcoord
                    data/measurement_ycoord
                    data/measurement_zcoord
                    data/measurement_time
                    data/measurement_values
                    params/reaction_rate'
  []
  [from_forward]
    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."}>>> = 'adjoint/inner_product data/objective_value'
    to_reporters<<<{"description": "List of the reporter names (object_name/value_name) to transfer the value to."}>>> = 'OptimizationReporter/grad_reaction_rate OptimizationReporter/objective_value'
  []
[]

[Reporters]
  [optInfo]
    type = OptimizationInfo<<<{"description": "Reports Optimization Output", "href": "../../../source/reporters/OptimizationInfo.html"}>>>
    items<<<{"description": "The information to output, if nothing is provided everything will be output."}>>> = 'current_iterate function_value gnorm'
  []
[]

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Optimize
  tao_solver = taobqnls
  petsc_options_iname = '-tao_gttol -tao_max_it'
  #petsc_options_value = '1e-5 100' #use this to get results for paper
  petsc_options_value = '1e-5 5'
  solve_on = 'NONE'
  verbose = true
[]
[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

Listing 2: Physics sub-application

[Mesh<<<{"href": "../../../syntax/Mesh/index.html"}>>>]
  [square]
    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"}>>> = 16
    ny<<<{"description": "Number of elements in the Y direction"}>>> = 16
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = 0
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 1
    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]
[]

[Reporters<<<{"href": "../../../syntax/Reporters/index.html"}>>>]
  [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."}>>> = 'reaction_rate'
    real_vector_values<<<{"description": "Values for vectors of reals."}>>> = '0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0' # Dummy
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = none
  []
  [data]
    type = OptimizationData<<<{"description": "Reporter to hold measurement and simulation data for optimization problems", "href": "../../../source/reporters/OptimizationData.html"}>>>
    variable<<<{"description": "Vector of variable names to sample at measurement points."}>>> = u
    objective_name<<<{"description": "Name of reporter value defining the objective."}>>> = objective_value
    measurement_file<<<{"description": "CSV file with measurement value and coordinates (value, x, y, z)."}>>> = forward_exact_csv_sample_0011.csv
    file_xcoord<<<{"description": "x coordinate column name from measurement_file csv being read in."}>>> = measurement_xcoord
    file_ycoord<<<{"description": "y coordinate column name from csv file being read in."}>>> = measurement_ycoord
    file_zcoord<<<{"description": "z coordinate column name from csv file being read in."}>>> = measurement_zcoord
    file_time<<<{"description": "time column name from csv file being read in."}>>> = measurement_time
    file_value<<<{"description": "measurement value column name from csv file being read in."}>>> = simulation_values
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = none
  []
[]

[Functions<<<{"href": "../../../syntax/Functions/index.html"}>>>]
  [rxn_func]
    type = ParameterMeshFunction<<<{"description": "Optimization function with parameters represented by a mesh and finite-element shape functions.", "href": "../../../source/functions/ParameterMeshFunction.html"}>>>
    exodus_mesh<<<{"description": "File containing parameter mesh."}>>> = parameter_mesh_out.e
    parameter_name<<<{"description": "Reporter or VectorPostprocessor vector containing parameter values."}>>> = params/reaction_rate
  []
[]

[Materials<<<{"href": "../../../syntax/Materials/index.html"}>>>]
  [ad_dc_prop]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../../source/materials/ParsedMaterial.html"}>>>
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '1 + u'
    coupled_variables<<<{"description": "Vector of variables used in the parsed function"}>>> = 'u'
    property_name<<<{"description": "Name of the parsed material property"}>>> = dc_prop
  []
  [ad_rxn_prop]
    type = ADGenericFunctionMaterial<<<{"description": "Material object for declaring properties that are populated by evaluation of Function object.", "href": "../../../source/materials/GenericFunctionMaterial.html"}>>>
    prop_values<<<{"description": "The corresponding names of the functions that are going to provide the values for the variables"}>>> = 'rxn_func'
    prop_names<<<{"description": "The names of the properties this material will have"}>>> = rxn_prop
  []
  #ADMatReaction includes a negative sign in residual evaluation, so we need to
  #reverse this with a negative reaction rate. However, we wanted the parameter
  #to remain positive, which is why there is one object to evaluate function
  #and another to flip it's sign for the kernel
  [ad_neg_rxn_prop]
    type = ADParsedMaterial<<<{"description": "Parsed expression Material.", "href": "../../../source/materials/ParsedMaterial.html"}>>>
    expression<<<{"description": "Parsed function (see FParser) expression for the parsed material"}>>> = '-rxn_prop'
    material_property_names<<<{"description": "Vector of material properties used in the parsed function"}>>> = 'rxn_prop'
    property_name<<<{"description": "Name of the parsed material property"}>>> = 'neg_rxn_prop'
  []
[]

[Kernels<<<{"href": "../../../syntax/Kernels/index.html"}>>>]
  [udot]
    type = ADTimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../../source/kernels/ADTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
  []
  [diff]
    type = ADMatDiffusion<<<{"description": "Diffusion equation kernel that takes an isotropic diffusivity from a material property", "href": "../../../source/kernels/ADMatDiffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    diffusivity<<<{"description": "The diffusivity value or material property"}>>> = dc_prop
  []
  [reaction]
    type = ADMatReaction<<<{"description": "Kernel representing the contribution of the PDE term $-L*v$, where $L$ is a reaction rate material property, $v$ is a scalar variable (nonlinear or coupled), and whose Jacobian contribution is calculated using automatic differentiation.", "href": "../../../source/kernels/ADMatReaction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    reaction_rate<<<{"description": "The reaction_rate used with the kernel."}>>> = neg_rxn_prop
  []
  [src]
    type = ADBodyForce<<<{"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
    value<<<{"description": "Coefficient to multiply by the body force term"}>>> = 1
  []
[]

[BCs<<<{"href": "../../../syntax/BCs/index.html"}>>>]
  [dirichlet]
    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 bottom'
    value<<<{"description": "Value of the BC"}>>> = 0
  []
[]

[Preconditioning<<<{"href": "../../../syntax/Preconditioning/index.html"}>>>]
  [nl0]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    nl_sys<<<{"description": "The nonlinear system whose linearization this preconditioner should be applied to."}>>> = 'nl0'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'lu'
  []
  [adjoint]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    nl_sys<<<{"description": "The nonlinear system whose linearization this preconditioner should be applied to."}>>> = 'adjoint'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-pc_type'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'lu'
  []
[]

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = TransientAndAdjoint
  forward_system = nl0
  adjoint_system = adjoint
  dt = 0.1
  end_time = 1

  nl_rel_tol = 1e-12
[]

[Problem<<<{"href": "../../../syntax/Problem/index.html"}>>>]
  nl_sys_names = 'nl0 adjoint'
  kernel_coverage_check = false
  skip_nl_system_check = true
[]

[Variables]
  [u_adjoint]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0
    solver_sys = adjoint
    outputs = none
  []
[]

[DiracKernels<<<{"href": "../../../syntax/DiracKernels/index.html"}>>>]
  [misfit]
    type = ReporterTimePointSource<<<{"description": "Apply a time dependent point load defined by Reporters.", "href": "../../../source/dirackernels/ReporterTimePointSource.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u_adjoint
    value_name<<<{"description": "reporter value name.  This uses the reporter syntax <reporter>/<name>."}>>> = data/misfit_values
    x_coord_name<<<{"description": "reporter x-coordinate name.  This uses the reporter syntax <reporter>/<name>."}>>> = data/measurement_xcoord
    y_coord_name<<<{"description": "reporter y-coordinate name.  This uses the reporter syntax <reporter>/<name>."}>>> = data/measurement_ycoord
    z_coord_name<<<{"description": "reporter z-coordinate name.  This uses the reporter syntax <reporter>/<name>."}>>> = data/measurement_zcoord
    time_name<<<{"description": "Name of vector-postprocessor or reporter vector containing time, default is assumed to be all 0s."}>>> = data/measurement_time
  []
[]

[VectorPostprocessors<<<{"href": "../../../syntax/VectorPostprocessors/index.html"}>>>]
  [adjoint]
    type = ElementOptimizationReactionFunctionInnerProduct<<<{"description": "Compute the gradient for reaction material inversion by taking the inner product of the forward and adjoint variables multiplied by the derivative of the optimization function with respect to the controllable parameters.", "href": "../../../source/vectorpostprocessors/ElementOptimizationReactionFunctionInnerProduct.html"}>>>
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = u_adjoint
    forward_variable<<<{"description": "Variable from the forward solution"}>>> = u
    function<<<{"description": "Optimization function."}>>> = rxn_func
    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
    outputs<<<{"description": "Vector of output names where you would like to restrict the output of variables(s) associated with this object"}>>> = none
  []
[]

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

[AuxKernels<<<{"href": "../../../syntax/AuxKernels/index.html"}>>>]
  [reaction_rate_aux]
    type = FunctionAux<<<{"description": "Auxiliary Kernel that creates and updates a field variable by sampling a function through space and time.", "href": "../../../source/auxkernels/FunctionAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = reaction_rate
    function<<<{"description": "The function to use as the value"}>>> = rxn_func
    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."}>>> = TIMESTEP_END
  []
[]

[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [u1]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = u
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0.25 0.25 0'
  []
  [u2]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = u
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '0.75 0.75 0'
  []
  [u3]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../source/postprocessors/PointValue.html"}>>>
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = u
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '1 1 0'
  []
[]

[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
  exodus<<<{"description": "Output the results using the default settings for Exodus output."}>>> = true
  console<<<{"description": "Output the results using the default settings for Console output"}>>> = false
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]

The top right figure in Figure 1 shows the rate found from the optimization process. Figure 2 shows how close the solution from the optimized reaction rate is from the synthetic measurement data. The diffusion-only initial guess is shown by the square data points.

Figure 1: Left: Exact reaction rate, simulation mesh, and measurement locations. Right: Optimized rection rate and parameter mesh.

Figure 2: Comparing simulated transient solution with exact, initial optimization guess, and optimized reaction rate at several measurement locations (measurement and optimized data are visually identical).

(modules/optimization/examples/diffusion_reaction/optimize.i)

[Optimization]
[]

[OptimizationReporter]
  type = ParameterMeshOptimization
  objective_name = objective_value
  parameter_names = 'reaction_rate'
  parameter_meshes = 'parameter_mesh_out.e'

  initial_condition = 0
  lower_bounds = 0

[]

[Reporters]
  [main]
    type = OptimizationData
    measurement_file = forward_exact_csv_sample_0011.csv
    file_xcoord = measurement_xcoord
    file_ycoord = measurement_ycoord
    file_zcoord = measurement_zcoord
    file_time = measurement_time
    file_value = simulation_values
  []
[]

[MultiApps]
  [forward]
    type = FullSolveMultiApp
    input_files = forward_and_adjoint.i
    execute_on = FORWARD
  []
[]

[Transfers]
  [to_forward]
    type = MultiAppReporterTransfer
    to_multi_app = forward
    from_reporters = 'main/measurement_xcoord
                      main/measurement_ycoord
                      main/measurement_zcoord
                      main/measurement_time
                      main/measurement_values
                      OptimizationReporter/reaction_rate'
    to_reporters = 'data/measurement_xcoord
                    data/measurement_ycoord
                    data/measurement_zcoord
                    data/measurement_time
                    data/measurement_values
                    params/reaction_rate'
  []
  [from_forward]
    type = MultiAppReporterTransfer
    from_multi_app = forward
    from_reporters = 'adjoint/inner_product data/objective_value'
    to_reporters = 'OptimizationReporter/grad_reaction_rate OptimizationReporter/objective_value'
  []
[]

[Reporters]
  [optInfo]
    type = OptimizationInfo
    items = 'current_iterate function_value gnorm'
  []
[]

[Executioner]
  type = Optimize
  tao_solver = taobqnls
  petsc_options_iname = '-tao_gttol -tao_max_it'
  #petsc_options_value = '1e-5 100' #use this to get results for paper
  petsc_options_value = '1e-5 5'
  solve_on = 'NONE'
  verbose = true
[]
[Outputs]
  csv = true
[]

(modules/optimization/examples/diffusion_reaction/forward_and_adjoint.i)

[Mesh]
  [square]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 16
    ny = 16
    xmin = 0
    xmax = 1
    ymin = 0
    ymax = 1
  []
[]

[Variables/u]
[]

[Reporters]
  [params]
    type = ConstantReporter
    real_vector_names = 'reaction_rate'
    real_vector_values = '0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0' # Dummy
    outputs = none
  []
  [data]
    type = OptimizationData
    variable = u
    objective_name = objective_value
    measurement_file = forward_exact_csv_sample_0011.csv
    file_xcoord = measurement_xcoord
    file_ycoord = measurement_ycoord
    file_zcoord = measurement_zcoord
    file_time = measurement_time
    file_value = simulation_values
    outputs = none
  []
[]

[Functions]
  [rxn_func]
    type = ParameterMeshFunction
    exodus_mesh = parameter_mesh_out.e
    parameter_name = params/reaction_rate
  []
[]

[Materials]
  [ad_dc_prop]
    type = ADParsedMaterial
    expression = '1 + u'
    coupled_variables = 'u'
    property_name = dc_prop
  []
  [ad_rxn_prop]
    type = ADGenericFunctionMaterial
    prop_values = 'rxn_func'
    prop_names = rxn_prop
  []
  #ADMatReaction includes a negative sign in residual evaluation, so we need to
  #reverse this with a negative reaction rate. However, we wanted the parameter
  #to remain positive, which is why there is one object to evaluate function
  #and another to flip it's sign for the kernel
  [ad_neg_rxn_prop]
    type = ADParsedMaterial
    expression = '-rxn_prop'
    material_property_names = 'rxn_prop'
    property_name = 'neg_rxn_prop'
  []
[]

[Kernels]
  [udot]
    type = ADTimeDerivative
    variable = u
  []
  [diff]
    type = ADMatDiffusion
    variable = u
    diffusivity = dc_prop
  []
  [reaction]
    type = ADMatReaction
    variable = u
    reaction_rate = neg_rxn_prop
  []
  [src]
    type = ADBodyForce
    variable = u
    value = 1
  []
[]

[BCs]
  [dirichlet]
    type = DirichletBC
    variable = u
    boundary = 'left bottom'
    value = 0
  []
[]

[Preconditioning]
  [nl0]
    type = SMP
    nl_sys = 'nl0'
    petsc_options_iname = '-pc_type'
    petsc_options_value = 'lu'
  []
  [adjoint]
    type = SMP
    nl_sys = 'adjoint'
    petsc_options_iname = '-pc_type'
    petsc_options_value = 'lu'
  []
[]

[Executioner]
  type = TransientAndAdjoint
  forward_system = nl0
  adjoint_system = adjoint
  dt = 0.1
  end_time = 1

  nl_rel_tol = 1e-12
[]

[Problem]
  nl_sys_names = 'nl0 adjoint'
  kernel_coverage_check = false
  skip_nl_system_check = true
[]

[Variables]
  [u_adjoint]
    initial_condition = 0
    solver_sys = adjoint
    outputs = none
  []
[]

[DiracKernels]
  [misfit]
    type = ReporterTimePointSource
    variable = u_adjoint
    value_name = data/misfit_values
    x_coord_name = data/measurement_xcoord
    y_coord_name = data/measurement_ycoord
    z_coord_name = data/measurement_zcoord
    time_name = data/measurement_time
  []
[]

[VectorPostprocessors]
  [adjoint]
    type = ElementOptimizationReactionFunctionInnerProduct
    variable = u_adjoint
    forward_variable = u
    function = rxn_func
    execute_on = ADJOINT_TIMESTEP_END
    outputs = none
  []
[]

[AuxVariables]
  [reaction_rate]
  []
[]

[AuxKernels]
  [reaction_rate_aux]
    type = FunctionAux
    variable = reaction_rate
    function = rxn_func
    execute_on = TIMESTEP_END
  []
[]

[Postprocessors]
  [u1]
    type = PointValue
    variable = u
    point = '0.25 0.25 0'
  []
  [u2]
    type = PointValue
    variable = u
    point = '0.75 0.75 0'
  []
  [u3]
    type = PointValue
    variable = u
    point = '1 1 0'
  []
[]

[Outputs]
  exodus = true
  console = false
  csv = true
[]