Merge pull request #28872 from maxnezdyur/line_obj

Gaussian Probe optimization

modules/optimization/doc/content/source/materials/MisfitReporterOffsetFunctionMaterial.md

@@ -0,0 +1,43 @@
+# MisfitReporterOffsetFunctionMaterial
+
+!syntax description /Materials/MisfitReporterOffsetFunctionMaterial
+
+## Overview
+
+The `MisfitReporterOffsetFunctionMaterial` material computes the misfit and its gradient with respect to a simulation variable , using a user-defined offset function to weight the contribution of measurement points. This material is particularly useful for inverse optimization problems where you want to match simulation results to experimental or observed data.  The material name given to the misfit is [!param](/Materials/MisfitReporterOffsetFunctionMaterial/property_name) and the misfit gradient material uses the naming convention [!param](/Materials/MisfitReporterOffsetFunctionMaterial/property_name)+`_gradient`.
+
+### Measurement data:
+
+Measurement values are provide by the reporter [!param](/Materials/MisfitReporterOffsetFunctionMaterial/measurement_value_name) and the corresponding measurement locations are given by a reporter containing either a vector of points [!param](/Materials/MisfitReporterOffsetFunctionMaterial/point_name) or three separate reporters each containing the x, y and z coordinates respectively of the points using [!param](/Materials/MisfitReporterOffsetFunctionMaterial/x_coord_name), [!param](/Materials/MisfitReporterOffsetFunctionMaterial/y_coord_name), and [!param](/Materials/MisfitReporterOffsetFunctionMaterial/z_coord_name).
+
+- Measurement values: $\mathbf{u}_m = \{u_{m,1}, u_{m,2}, \ldots, u_{m,n}\}$
+- Measurement locations: $\mathbf{p}_m = \{\mathbf{p}_{m,1}, \mathbf{p}_{m,2}, \ldots, \mathbf{p}_{m,n}\}$, where $\mathbf{p}_{m,i} = (x_{m,i}, y_{m,i}, z_{m,i})$ for 3D data
+- Simulation variable: $u(\mathbf{x})$
+- Offset function: $g(\mathbf{p}, \mathbf{x}; \mathbf{\theta})$, where $\mathbf{\theta}$ are the parameters defining the offset function
+- Simulation point: $\mathbf{x}$
+
+### Misfit and gradient calculation:
+
+The `misfit` at a given simulation point $\mathbf{x}$ is calculated as:
+
+$m(\mathbf{x}) = \sum_{i=1}^{n} \left( u_{m,i} g(\mathbf{p}_{m,i}, \mathbf{x}; \mathbf{\theta}) - u(\mathbf{x}) g(\mathbf{p}_{m,i}, \mathbf{x}; \mathbf{\theta}) \right)^2,$
+
+and the `misfit_gradient` with respect to the simulation variable $u(\mathbf{x})$ is:
+
+$\frac{\partial m(\mathbf{x})}{\partial u(\mathbf{x})} = -2 \sum_{i=1}^{n} g(\mathbf{p}_{m,i}, \mathbf{x}; \mathbf{\theta}) \left( u_{m,i} g(\mathbf{p}_{m,i}, \mathbf{x}; \mathbf{\theta}) - u(\mathbf{x}) g(\mathbf{p}_{m,i}, \mathbf{x}; \mathbf{\theta}) \right).$
+
+These equations represent the misfit value and its gradient computed by the `MisfitReporterOffsetFunctionMaterial` material at each quadrature point in the simulation. The misfit value quantifies the discrepancy between the measured data and the simulation variable, while the misfit gradient provides the sensitivity of the misfit with respect to changes in the simulation variable.
+
+## Example Input File Syntax
+
+The following input file contains an example of how to use this material to define gaussian shaped measurement points in a heat source inversion problem where the misfit is used in the computation of the objective function and the misfit_gradient is used to compute the right hand side in the adjoint problem.
+
+!listing test/tests/optimizationreporter/function_misfit/forward_and_adjoint.i block=Materials/beam
+
+## Syntax parameters and inputs
+
+!syntax parameters /Materials/MisfitReporterOffsetFunctionMaterial
+
+!syntax inputs /Materials/MisfitReporterOffsetFunctionMaterial
+
+!syntax children /Materials/MisfitReporterOffsetFunctionMaterial

modules/optimization/doc/content/source/materials/ReporterOffsetFunctionMaterial.md

@@ -0,0 +1,25 @@
+# ReporterOffsetFunctionMaterial
+
+!syntax description /Materials/ReporterOffsetFunctionMaterial
+
+This creates a material that is the sum of a function shifted by a set of points contained in a reporter.  The reporter can contain points by specifying [!param](/Materials/ReporterOffsetFunctionMaterial/point_name) or by three seperate reporters containing the (x,y,z) coordinates of the points using [!param](/Materials/ReporterOffsetFunctionMaterial/x_coord_name), [!param](/Materials/ReporterOffsetFunctionMaterial/y_coord_name), and [!param](/Materials/ReporterOffsetFunctionMaterial/z_coord_name). This can be useful for creating a field containing multiple sources, as shown in [offset] and described in the equation below.
+
+
+!equation
+\tilde f(\mathbf{x}) = \sum_{i=1}^{n} \left(f(\mathbf{x} - \mathbf{p}_{offset}) \right)
+
+!media large_media/framework/materials/offset_func.png
+       id=offset
+       caption=Field with multiple offsets of (left) gaussian function and (right) constant circle function.
+
+## Example Input File Syntax
+
+In this example, `ReporterOffsetFunctionMaterial` is used to define the two materials shown in [offset]. The value at a point is the sum of all the shifted functions.  In this example, the gaussian and constant circle functions are defined at (x,y,z)=(0,0,0), shown by the source in the lower left corner of each domain.  This function is then shifted according to the offset vectors given by the Reporter.
+
+!listing test/tests/optimizationreporter/reporter_offset/reporter_offset_mat.i block=Materials/multiple_sources
+
+!syntax parameters /Materials/ReporterOffsetFunctionMaterial
+
+!syntax inputs /Materials/ReporterOffsetFunctionMaterial
+
+!syntax children /Materials/ReporterOffsetFunctionMaterial

modules/optimization/include/materials/MisfitReporterOffsetFunctionMaterial.h

@@ -0,0 +1,54 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "ReporterOffsetFunctionMaterial.h"
+
+template <bool>
+class MisfitReporterOffsetFunctionMaterialTempl;
+typedef MisfitReporterOffsetFunctionMaterialTempl<false> MisfitReporterOffsetFunctionMaterial;
+typedef MisfitReporterOffsetFunctionMaterialTempl<true> ADMisfitReporterOffsetFunctionMaterial;
+
+/**
+ * This class creates a misfit and misfit gradient material that can be used for
+ * optimizing measurements weighted by offset functions.
+ */
+template <bool is_ad>
+class MisfitReporterOffsetFunctionMaterialTempl : public ReporterOffsetFunctionMaterialTempl<is_ad>
+{
+public:
+  static InputParameters validParams();
+
+  MisfitReporterOffsetFunctionMaterialTempl<is_ad>(const InputParameters & parameters);
+
+protected:
+  virtual void computeQpProperties() override;
+
+  /// Simulation variable
+  const GenericVariableValue<is_ad> & _sim_var;
+
+  /// Gradient of misfit with respect to the simulation variable
+  GenericMaterialProperty<Real, is_ad> & _mat_prop_gradient;
+
+  /// values at each xyz coordinate
+  const std::vector<Real> & _measurement_values;
+
+  Real computeWeighting(const Point & point_shift);
+
+  using ReporterOffsetFunctionMaterialTempl<is_ad>::_prop_name;
+  using ReporterOffsetFunctionMaterialTempl<is_ad>::_qp;
+  using ReporterOffsetFunctionMaterialTempl<is_ad>::_material;
+  using ReporterOffsetFunctionMaterialTempl<is_ad>::_points;
+  using ReporterOffsetFunctionMaterialTempl<is_ad>::_read_in_points;
+  using ReporterOffsetFunctionMaterialTempl<is_ad>::_coordx;
+  using ReporterOffsetFunctionMaterialTempl<is_ad>::_coordy;
+  using ReporterOffsetFunctionMaterialTempl<is_ad>::_coordz;
+  using ReporterOffsetFunctionMaterialTempl<is_ad>::computeOffsetFunction;
+};

modules/optimization/include/materials/ReporterOffsetFunctionMaterial.h

@@ -0,0 +1,72 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#pragma once
+
+#include "Material.h"
+#include "Function.h"
+#include "ReporterInterface.h"
+
+template <bool>
+class ReporterOffsetFunctionMaterialTempl;
+typedef ReporterOffsetFunctionMaterialTempl<false> ReporterOffsetFunctionMaterial;
+typedef ReporterOffsetFunctionMaterialTempl<true> ADReporterOffsetFunctionMaterial;
+
+/**
+ * This class defines a material with an associated offset function.
+ * It is designed to be used in conjunction with a reporter, allowing
+ * for flexible and customizable material properties in a simulation.
+ */
+template <bool is_ad>
+class ReporterOffsetFunctionMaterialTempl : public Material, public ReporterInterface
+{
+public:
+  static InputParameters validParams();
+
+  ReporterOffsetFunctionMaterialTempl(const InputParameters & parameters);
+
+protected:
+  virtual void computeQpProperties() override;
+
+  /// Base name of the material system
+  const std::string _prop_name;
+
+  /// Misfit value at each quadrature point
+  GenericMaterialProperty<Real, is_ad> & _material;
+
+  /// bool if data format read in is points
+  const bool _read_in_points;
+  /**
+   * Reporter offset locations for function
+   */
+  // @{
+  /// x coordinate
+  const std::vector<Real> & _coordx;
+  /// y coordinate
+  const std::vector<Real> & _coordy;
+  ///z coordinate
+  const std::vector<Real> & _coordz;
+  ///xyz point
+  const std::vector<Point> & _points;
+  // @}
+
+  /// The function being used for evaluation
+  const Function & _func;
+  /**
+   * Calculates the value of the function at the given offset point.
+   * @param point_offset The offset point to shift the function evaluation by.
+   * @return The value of the function at the shifted location.
+   */
+  Real computeOffsetFunction(const Point & point_offset);
+
+private:
+  /// convenience vectors (these are not const because reporters can change their size)
+  std::vector<Real> _zeros_vec;
+  std::vector<Point> _zeros_pts;
+};

modules/optimization/src/materials/MisfitReporterOffsetFunctionMaterial.C

@@ -0,0 +1,78 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "MisfitReporterOffsetFunctionMaterial.h"
+
+#include "libmesh/libmesh_common.h"
+
+registerMooseObject("OptimizationApp", MisfitReporterOffsetFunctionMaterial);
+registerMooseObject("OptimizationApp", ADMisfitReporterOffsetFunctionMaterial);
+
+template <bool is_ad>
+InputParameters
+MisfitReporterOffsetFunctionMaterialTempl<is_ad>::validParams()
+{
+  InputParameters params = ReporterOffsetFunctionMaterial::validParams();
+  params.addClassDescription(
+      "Computes the misfit and misfit gradient materials for inverse optimizations problems.");
+
+  params.addRequiredCoupledVar("forward_variable",
+                               "Variable that is being used for the forward simulation.");
+  params.addRequiredParam<ReporterName>("measurement_value_name",
+                                        "Reporter with measurement data.");
+  return params;
+}
+
+template <bool is_ad>
+MisfitReporterOffsetFunctionMaterialTempl<is_ad>::MisfitReporterOffsetFunctionMaterialTempl(
+    const InputParameters & parameters)
+  : ReporterOffsetFunctionMaterialTempl<is_ad>(parameters),
+    _sim_var(this->template coupledGenericValue<is_ad>("forward_variable")),
+    _mat_prop_gradient(
+        this->template declareGenericProperty<Real, is_ad>(_prop_name + "_gradient")),
+    _measurement_values(this->template getReporterValue<std::vector<Real>>(
+        "measurement_value_name", REPORTER_MODE_REPLICATED))
+{
+}
+
+template <bool is_ad>
+void
+MisfitReporterOffsetFunctionMaterialTempl<is_ad>::computeQpProperties()
+{
+  _material[_qp] = 0.0;
+  _mat_prop_gradient[_qp] = 0.0;
+  auto num_pts = _read_in_points ? _points.size() : _coordx.size();
+  if (!_read_in_points)
+    mooseAssert((_coordx.size() == _coordy.size()) && (_coordx.size() == _coordz.size()),
+                "Size of the coordinate offsets don't match.");
+
+  mooseAssert(num_pts == _measurement_values.size(),
+              "Number of offsets doesn't match the number of measurements.");
+
+  for (const auto idx : make_range(num_pts))
+  {
+    const Point offset =
+        _read_in_points ? _points[idx] : Point(_coordx[idx], _coordy[idx], _coordz[idx]);
+
+    const Real measurement_value = _measurement_values[idx];
+    const auto simulation_value = _sim_var[_qp];
+
+    // Compute weighting function
+    const Real weighting = computeOffsetFunction(offset);
+
+    // Computed weighted misfit and gradient materials
+    _material[_qp] +=
+        Utility::pow<2>(weighting) * Utility::pow<2>(measurement_value - simulation_value);
+    _mat_prop_gradient[_qp] -=
+        2.0 * Utility::pow<2>(weighting) * (measurement_value - simulation_value);
+  }
+}
+
+template class MisfitReporterOffsetFunctionMaterialTempl<true>;
+template class MisfitReporterOffsetFunctionMaterialTempl<false>;

modules/optimization/src/materials/ReporterOffsetFunctionMaterial.C

@@ -0,0 +1,89 @@
+//* This file is part of the MOOSE framework
+//* https://www.mooseframework.org
+//*
+//* All rights reserved, see COPYRIGHT for full restrictions
+//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
+//*
+//* Licensed under LGPL 2.1, please see LICENSE for details
+//* https://www.gnu.org/licenses/lgpl-2.1.html
+
+#include "ReporterOffsetFunctionMaterial.h"
+
+registerMooseObject("OptimizationApp", ReporterOffsetFunctionMaterial);
+registerMooseObject("OptimizationApp", ADReporterOffsetFunctionMaterial);
+
+template <bool is_ad>
+InputParameters
+ReporterOffsetFunctionMaterialTempl<is_ad>::validParams()
+{
+  InputParameters params = Material::validParams();
+  params.addClassDescription("Compute the sum of a function offset by a set of points.");
+  params.addRequiredParam<FunctionName>("function", "The weighting function.");
+
+  params.addParam<ReporterName>("x_coord_name", "Reporter with x-coordinate data.");
+  params.addParam<ReporterName>("y_coord_name", "Reporter with y-coordinate data.");
+  params.addParam<ReporterName>("z_coord_name", "Reporter with z-coordinate data.");
+  params.addParam<ReporterName>("point_name",
+                                "Reporter with point data. "
+                                "<reporter>/<name>.");
+  params.addRequiredParam<std::string>("property_name", "Material property base name");
+  params.addParamNamesToGroup("point_name x_coord_name y_coord_name z_coord_name",
+                              "Offset locations for function evaluations");
+  return params;
+}
+
+template <bool is_ad>
+ReporterOffsetFunctionMaterialTempl<is_ad>::ReporterOffsetFunctionMaterialTempl(
+    const InputParameters & parameters)
+  : Material(parameters),
+    ReporterInterface(this),
+    _prop_name(getParam<std::string>("property_name")),
+    _material(declareGenericProperty<Real, is_ad>(_prop_name)),
+    _read_in_points(isParamValid("point_name")),
+    _coordx(isParamValid("x_coord_name")
+                ? getReporterValue<std::vector<Real>>("x_coord_name", REPORTER_MODE_REPLICATED)
+                : _zeros_vec),
+    _coordy(isParamValid("y_coord_name")
+                ? getReporterValue<std::vector<Real>>("y_coord_name", REPORTER_MODE_REPLICATED)
+                : _zeros_vec),
+    _coordz(isParamValid("z_coord_name")
+                ? getReporterValue<std::vector<Real>>("z_coord_name", REPORTER_MODE_REPLICATED)
+                : _zeros_vec),
+    _points(_read_in_points
+                ? getReporterValue<std::vector<Point>>("point_name", REPORTER_MODE_REPLICATED)
+                : _zeros_pts),
+    _func(getFunction("function"))
+{
+  if (isParamValid("point_name") == (isParamValid("x_coord_name") && isParamValid("y_coord_name") &&
+                                     isParamValid("z_coord_name")))
+    paramError("Either supply x,y, and z reporters or a point reporter.");
+}
+
+template <bool is_ad>
+void
+ReporterOffsetFunctionMaterialTempl<is_ad>::computeQpProperties()
+{
+  _material[_qp] = 0;
+  auto num_pts = _read_in_points ? _points.size() : _coordx.size();
+  if (!_read_in_points)
+    mooseAssert((_coordx.size() == _coordy.size()) && (_coordx.size() == _coordz.size()),
+                "Size of the coordinate offsets don't match.");
+
+  for (const auto idx : make_range(num_pts))
+  {
+
+    Point offset = _read_in_points ? _points[idx] : Point(_coordx[idx], _coordy[idx], _coordz[idx]);
+
+    _material[_qp] += computeOffsetFunction(offset);
+  }
+}
+
+template <bool is_ad>
+Real
+ReporterOffsetFunctionMaterialTempl<is_ad>::computeOffsetFunction(const Point & point_offset)
+{
+  return _func.value(_t, _q_point[_qp] - point_offset);
+}
+
+template class ReporterOffsetFunctionMaterialTempl<true>;
+template class ReporterOffsetFunctionMaterialTempl<false>;