functional_objective.hpp

Go to the documentation of this file.

// Copyright Lawrence Livermore National Security, LLC and
// other Smith Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include "smith/physics/scalar_objective.hpp"
#include "smith/physics/mesh.hpp"
#include "smith/numerics/functional/shape_aware_functional.hpp"
#include "smith/physics/state/finite_element_state.hpp"
#include "smith/physics/state/finite_element_dual.hpp"
 
namespace smith {
 
template <int spatial_dim, typename parameters = Parameters<>,
          typename parameter_indices = std::make_integer_sequence<int, parameters::n>>
class FunctionalObjective;
 
template <int spatial_dim, typename... InputSpaces, int... parameter_indices>
class FunctionalObjective<spatial_dim, Parameters<InputSpaces...>, std::integer_sequence<int, parameter_indices...>>
    : public ScalarObjective {
 public:
  using SpacesT = std::vector<const mfem::ParFiniteElementSpace*>;  
 
  using ShapeDispSpace = H1<1, spatial_dim>;  
 
  FunctionalObjective(const std::string& physics_name, std::shared_ptr<Mesh> mesh, const SpacesT& input_mfem_spaces)
      : ScalarObjective(physics_name), mesh_(mesh)
  {
    std::array<const mfem::ParFiniteElementSpace*, sizeof...(InputSpaces)> mfem_spaces;
 
    SLIC_ERROR_ROOT_IF(
        sizeof...(InputSpaces) != input_mfem_spaces.size(),
        std::format("{} parameter spaces given in the template argument but {} parameter names were supplied.",
                    sizeof...(InputSpaces), input_mfem_spaces.size()));
 
    if constexpr (sizeof...(InputSpaces) > 0) {
      for_constexpr<sizeof...(InputSpaces)>([&](auto i) { mfem_spaces[i] = input_mfem_spaces[i]; });
    }
 
    const auto& shape_disp_space = mesh_->shapeDisplacementSpace();
 
    objective_ =
        std::make_unique<ShapeAwareFunctional<ShapeDispSpace, double(InputSpaces...)>>(&shape_disp_space, mfem_spaces);
  }
 
  template <typename FuncOfTimeSpaceAndParams, int... all_params>
  void addBodyIntegralImpl(std::string body_name, const FuncOfTimeSpaceAndParams& qfunction,
                           std::integer_sequence<int, all_params...>)
  {
    const double* dt = &dt_;
    const size_t* cycle = &cycle_;
    objective_->AddDomainIntegral(
        smith::Dimension<spatial_dim>{}, smith::DependsOn<all_params...>{},
        [dt, cycle, qfunction](double time, auto X, auto... params) {
          return qfunction(TimeInfo(time, *dt, *cycle), X, params...);
        },
        mesh_->domain(body_name));
  }
 
  template <int... active_parameters, typename FuncOfTimeSpaceAndParams>
  void addBodyIntegral(DependsOn<active_parameters...>, std::string body_name,
                       const FuncOfTimeSpaceAndParams& qfunction)
  {
    addBodyIntegralImpl(body_name, qfunction, std::integer_sequence<int, active_parameters...>{});
  }
 
  template <typename FuncOfTimeSpaceAndParams>
  void addBodyIntegral(std::string body_name, const FuncOfTimeSpaceAndParams& qfunction)
  {
    addBodyIntegralImpl(body_name, qfunction, std::make_integer_sequence<int, sizeof...(InputSpaces)>{});
  }
 
  template <typename FuncOfTimeSpaceAndParams, int... all_params>
  void addBoundaryIntegralImpl(std::string boundary_name, const FuncOfTimeSpaceAndParams& qfunction,
                               std::integer_sequence<int, all_params...>)
  {
    const double* dt = &dt_;
    const size_t* cycle = &cycle_;
    objective_->AddBoundaryIntegral(
        smith::Dimension<spatial_dim>{}, smith::DependsOn<all_params...>{},
        [dt, cycle, qfunction](double time, auto X, auto... params) {
          return qfunction(TimeInfo(time, *dt, *cycle), X, params...);
        },
        mesh_->domain(boundary_name));
  }
 
  template <int... active_parameters, typename FuncOfTimeSpaceAndParams>
  void addBoundaryIntegral(DependsOn<active_parameters...>, std::string boundary_name,
                           const FuncOfTimeSpaceAndParams& qfunction)
  {
    addBoundaryIntegralImpl(boundary_name, qfunction, std::integer_sequence<int, active_parameters...>{});
  }
 
  template <typename FuncOfTimeSpaceAndParams>
  void addBoundaryIntegral(std::string boundary_name, const FuncOfTimeSpaceAndParams& qfunction)
  {
    addBoundaryIntegralImpl(boundary_name, qfunction, std::make_integer_sequence<int, sizeof...(InputSpaces)>{});
  }
 
  virtual double evaluate(TimeInfo time_info, ConstFieldPtr shape_disp,
                          const std::vector<ConstFieldPtr>& fields) const override
  {
    dt_ = time_info.dt();
    cycle_ = time_info.cycle();
 
    return evaluateObjective(std::make_integer_sequence<int, sizeof...(parameter_indices)>{}, time_info.time(),
                             shape_disp, fields);
  }
 
  virtual mfem::Vector gradient(TimeInfo time_info, ConstFieldPtr shape_disp, const std::vector<ConstFieldPtr>& fields,
                                size_t field_ordinal) const override
  {
    dt_ = time_info.dt();
    cycle_ = time_info.cycle();
 
    auto grads = gradientEvaluators(std::make_integer_sequence<int, sizeof...(parameter_indices)>{}, time_info.time(),
                                    shape_disp, fields);
    auto g = smith::get<DERIVATIVE>(grads[field_ordinal](time_info.time(), shape_disp, fields));
    return *assemble(g);
  }
 
  virtual mfem::Vector mesh_coordinate_gradient(TimeInfo time_info, ConstFieldPtr shape_disp,
                                                const std::vector<ConstFieldPtr>& fields) const override
  {
    dt_ = time_info.dt();
    cycle_ = time_info.cycle();
 
    auto g = smith::get<DERIVATIVE>(
        (*objective_)(DifferentiateWRT<0>{}, time_info.time(), *shape_disp, *fields[parameter_indices]...));
    return *assemble(g);
  }
 
 private:
  template <int... i>
  auto evaluateObjective(std::integer_sequence<int, i...>, double time, ConstFieldPtr shape_disp,
                         const std::vector<ConstFieldPtr>& fs) const
  {
    return (*objective_)(time, *shape_disp, *fs[i]...);
  }
 
  template <int... i>
  auto gradientEvaluators(std::integer_sequence<int, i...>, double time, ConstFieldPtr shape_disp,
                          const std::vector<ConstFieldPtr>& fs) const
  {
    using JacFuncType = std::function<decltype((*objective_)(DifferentiateWRT<1>{}, time, *shape_disp, *fs[i]...))(
        double, ConstFieldPtr, const std::vector<ConstFieldPtr>&)>;
    return std::array<JacFuncType, sizeof...(i)>{
        [this](double _time, ConstFieldPtr _shape_disp, const std::vector<ConstFieldPtr>& _fs) {
          return (*objective_)(DifferentiateWRT<i + 1>{}, _time, *_shape_disp, *_fs[i]...);
        }...};
  }
 
  mutable double dt_ = std::numeric_limits<double>::max();
 
  mutable size_t cycle_ = 0;
 
  std::shared_ptr<Mesh> mesh_;
 
  std::unique_ptr<ShapeAwareFunctional<ShapeDispSpace, double(InputSpaces...)>> objective_;
};
 
}  // namespace smith