time_discretized_weak_form.hpp

Go to the documentation of this file.

// Copyright (c) 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/functional_weak_form.hpp"
#include "smith/physics/mesh.hpp"
#include "smith/differentiable_numerics/field_state.hpp"
#include "smith/differentiable_numerics/time_integration_rule.hpp"
 
namespace smith {
 
template <int spatial_dim, typename OutputSpace, typename inputs = Parameters<>>
class TimeDiscretizedWeakForm;
 
template <int spatial_dim, typename OutputSpace, typename... InputSpaces>
class TimeDiscretizedWeakForm<spatial_dim, OutputSpace, Parameters<InputSpaces...>>
    : public FunctionalWeakForm<spatial_dim, OutputSpace, Parameters<InputSpaces...>> {
 public:
  using Base = FunctionalWeakForm<spatial_dim, OutputSpace, Parameters<InputSpaces...>>;  
 
  TimeDiscretizedWeakForm(std::string physics_name, std::shared_ptr<Mesh> mesh,
                          const mfem::ParFiniteElementSpace& output_mfem_space,
                          const typename Base::SpacesT& input_mfem_spaces)
      : Base(physics_name, mesh, output_mfem_space, input_mfem_spaces)
  {
  }
 
  template <int... active_parameters, typename BodyIntegralType>
  void addBodyIntegral(DependsOn<active_parameters...> depends_on, std::string body_name, BodyIntegralType integrand)
  {
    const double* dt = &this->dt_;
    const size_t* cycle = &this->cycle_;
    Base::addBodyIntegral(depends_on, body_name, [dt, cycle, integrand](double t, auto X, auto... inputs) {
      TimeInfo time_info(t, *dt, *cycle);
      return integrand(time_info, X, inputs...);
    });
  }
 
  template <typename BodyIntegralType>
  void addBodyIntegral(std::string body_name, BodyIntegralType integrand)
  {
    constexpr int num_inputs = sizeof...(InputSpaces);
    addBodyIntegralWithAllParams(body_name, integrand, std::make_integer_sequence<int, num_inputs>{});
  }
 
  template <int... active_parameters, typename BodyLoadType>
  void addBodySource(DependsOn<active_parameters...> depends_on, std::string body_name, BodyLoadType load_function)
  {
    addBodyIntegral(depends_on, body_name, [load_function](auto t_info, auto X, auto... inputs) {
      return smith::tuple{-load_function(t_info, get<VALUE>(X), get<VALUE>(inputs)...), smith::zero{}};
    });
  }
 
  template <typename BodyLoadType>
  void addBodySource(std::string body_name, BodyLoadType load_function)
  {
    constexpr int num_inputs = sizeof...(InputSpaces);
    addBodySourceWithAllParams(body_name, load_function, std::make_integer_sequence<int, num_inputs>{});
  }
 
  template <int... active_parameters, typename BoundaryIntegralType>
  void addBoundaryIntegral(DependsOn<active_parameters...> depends_on, std::string boundary_name,
                           BoundaryIntegralType integrand)
  {
    const double* dt = &this->dt_;
    const size_t* cycle = &this->cycle_;
    Base::addBoundaryIntegral(depends_on, boundary_name, [dt, cycle, integrand](double t, auto X, auto... inputs) {
      TimeInfo time_info(t, *dt, *cycle);
      return integrand(time_info, X, inputs...);
    });
  }
 
  template <typename BoundaryIntegralType>
  void addBoundaryIntegral(std::string boundary_name, BoundaryIntegralType integrand)
  {
    constexpr int num_inputs = sizeof...(InputSpaces);
    addBoundaryIntegralWithAllParams(boundary_name, integrand, std::make_integer_sequence<int, num_inputs>{});
  }
 
  template <int... active_parameters, typename BoundaryFluxType>
  void addBoundaryFlux(DependsOn<active_parameters...> depends_on, std::string boundary_name,
                       BoundaryFluxType flux_function)
  {
    const double* dt = &this->dt_;
    const size_t* cycle = &this->cycle_;
    Base::addBoundaryFlux(depends_on, boundary_name,
                          [dt, cycle, flux_function](double t, auto X, auto n, auto... inputs) {
                            TimeInfo time_info(t, *dt, *cycle);
                            return flux_function(time_info, X, n, inputs...);
                          });
  }
 
  template <typename BoundaryFluxType>
  void addBoundaryFlux(std::string boundary_name, BoundaryFluxType flux_function)
  {
    constexpr int num_inputs = sizeof...(InputSpaces);
    addBoundaryFluxWithAllParams(boundary_name, flux_function, std::make_integer_sequence<int, num_inputs>{});
  }
 
  template <int... active_parameters, typename InteriorIntegralType>
  void addInteriorBoundaryIntegral(DependsOn<active_parameters...> depends_on, std::string interior_name,
                                   InteriorIntegralType integrand)
  {
    const double* dt = &this->dt_;
    const size_t* cycle = &this->cycle_;
    Base::addInteriorBoundaryIntegral(depends_on, interior_name,
                                      [dt, cycle, integrand](double t, auto X, auto... inputs) {
                                        TimeInfo time_info(t, *dt, *cycle);
                                        return integrand(time_info, X, inputs...);
                                      });
  }
 
  template <typename InteriorIntegralType>
  void addInteriorBoundaryIntegral(std::string interior_name, InteriorIntegralType integrand)
  {
    constexpr int num_inputs = sizeof...(InputSpaces);
    addInteriorBoundaryIntegralWithAllParams(interior_name, integrand, std::make_integer_sequence<int, num_inputs>{});
  }
 
 private:
  template <typename BodyIntegralType, int... all_params>
  void addBodyIntegralWithAllParams(std::string body_name, BodyIntegralType integrand,
                                    std::integer_sequence<int, all_params...>)
  {
    addBodyIntegral(DependsOn<all_params...>{}, body_name, integrand);
  }
 
  template <typename BodyLoadType, int... all_params>
  void addBodySourceWithAllParams(std::string body_name, BodyLoadType load_function,
                                  std::integer_sequence<int, all_params...>)
  {
    addBodySource(DependsOn<all_params...>{}, body_name, load_function);
  }
 
  template <typename BoundaryIntegralType, int... all_params>
  void addBoundaryIntegralWithAllParams(std::string boundary_name, BoundaryIntegralType integrand,
                                        std::integer_sequence<int, all_params...>)
  {
    addBoundaryIntegral(DependsOn<all_params...>{}, boundary_name, integrand);
  }
 
  template <typename InteriorIntegralType, int... all_params>
  void addInteriorBoundaryIntegralWithAllParams(std::string interior_name, InteriorIntegralType integrand,
                                                std::integer_sequence<int, all_params...>)
  {
    addInteriorBoundaryIntegral(DependsOn<all_params...>{}, interior_name, integrand);
  }
 
  template <typename BoundaryFluxType, int... all_params>
  void addBoundaryFluxWithAllParams(std::string boundary_name, BoundaryFluxType flux_function,
                                    std::integer_sequence<int, all_params...>)
  {
    addBoundaryFlux(DependsOn<all_params...>{}, boundary_name, flux_function);
  }
};
 
class SecondOrderTimeDiscretizedWeakForms {
 public:
  std::shared_ptr<WeakForm> time_discretized_weak_form;  
  std::shared_ptr<WeakForm> final_reaction_weak_form;    
};
 
template <int spatial_dim, typename OutputSpace, typename inputs = Parameters<>>
class SecondOrderTimeDiscretizedWeakForm;
 
template <int spatial_dim, typename OutputSpace, typename TrialInputSpace, typename... InputSpaces>
class SecondOrderTimeDiscretizedWeakForm<spatial_dim, OutputSpace, Parameters<TrialInputSpace, InputSpaces...>>
    : public SecondOrderTimeDiscretizedWeakForms {
 public:
  static constexpr int NUM_STATE_VARS = 4;  
 
  using TimeDiscretizedWeakFormT =
      TimeDiscretizedWeakForm<spatial_dim, OutputSpace, Parameters<TrialInputSpace, InputSpaces...>>;
  using FinalReactionFormT = TimeDiscretizedWeakForm<spatial_dim, OutputSpace, Parameters<InputSpaces...>>;
 
  SecondOrderTimeDiscretizedWeakForm(std::string physics_name, std::shared_ptr<Mesh> mesh,
                                     ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule,
                                     const mfem::ParFiniteElementSpace& output_mfem_space,
                                     const typename TimeDiscretizedWeakFormT::SpacesT& input_mfem_spaces)
      : time_rule_(time_rule)
  {
    time_discretized_weak_form_ =
        std::make_shared<TimeDiscretizedWeakFormT>(physics_name, mesh, output_mfem_space, input_mfem_spaces);
    time_discretized_weak_form = time_discretized_weak_form_;
 
    typename TimeDiscretizedWeakFormT::SpacesT trial_removed_spaces(std::next(input_mfem_spaces.begin()),
                                                                    input_mfem_spaces.end());
    final_reaction_weak_form_ =
        std::make_shared<FinalReactionFormT>(physics_name, mesh, output_mfem_space, trial_removed_spaces);
    final_reaction_weak_form = final_reaction_weak_form_;
  }
 
  template <int... active_parameters, typename BodyIntegralType>
  void addBodyIntegral(DependsOn<active_parameters...> /*depends_on*/, std::string body_name,
                       BodyIntegralType integrand)
  {
    auto time_rule = time_rule_;
    time_discretized_weak_form_->addBodyIntegral(
        DependsOn<0, 1, 2, 3, NUM_STATE_VARS + active_parameters...>{}, body_name,
        [integrand, time_rule](const TimeInfo& t, auto X, auto U, auto U_old, auto U_dot_old, auto U_dot_dot_old,
                               auto... inputs) {
          return integrand(t, X, time_rule.value(t, U, U_old, U_dot_old, U_dot_dot_old),
                           time_rule.dot(t, U, U_old, U_dot_old, U_dot_dot_old),
                           time_rule.ddot(t, U, U_old, U_dot_old, U_dot_dot_old), inputs...);
        });
    final_reaction_weak_form_->addBodyIntegral(DependsOn<0, 1, 2, NUM_STATE_VARS - 1 + active_parameters...>{},
                                               body_name, integrand);
  }
 
  template <typename BodyIntegralType>
  void addBodyIntegral(std::string body_name, BodyIntegralType integrand)
  {
    addBodyIntegral(DependsOn<>{}, body_name, integrand);
  }
 
  template <int... active_parameters, typename BodyLoadType>
  void addBodySource(DependsOn<active_parameters...> /*depends_on*/, std::string body_name, BodyLoadType load_function)
  {
    auto time_rule = time_rule_;
    time_discretized_weak_form_->addBodyIntegral(
        DependsOn<0, 1, 2, 3, NUM_STATE_VARS + active_parameters...>{}, body_name,
        [load_function, time_rule](const TimeInfo& t, auto X, auto U, auto U_old, auto U_dot_old, auto U_dot_dot_old,
                                   auto... inputs) {
          return smith::tuple{
              -load_function(t.time(), get<VALUE>(X),
                             get<VALUE>(time_rule.value(t, U, U_old, U_dot_old, U_dot_dot_old)),
                             get<VALUE>(time_rule.dot(t, U, U_old, U_dot_old, U_dot_dot_old)),
                             get<VALUE>(time_rule.ddot(t, U, U_old, U_dot_old, U_dot_dot_old)), get<VALUE>(inputs)...),
              smith::zero{}};
        });
    final_reaction_weak_form_->addBodyIntegral(
        DependsOn<0, 1, 2, NUM_STATE_VARS - 1 + active_parameters...>{}, body_name,
        [load_function](const TimeInfo& t, auto X, auto... inputs) {
          return smith::tuple{-load_function(t.time(), get<VALUE>(X), get<VALUE>(inputs)...), smith::zero{}};
        });
  }
 
  template <typename BodyLoadType>
  void addBodySource(std::string body_name, BodyLoadType load_function)
  {
    addBodySource(DependsOn<>{}, body_name, load_function);
  }
 
 private:
  std::shared_ptr<TimeDiscretizedWeakFormT> time_discretized_weak_form_;
  std::shared_ptr<FinalReactionFormT> final_reaction_weak_form_;
  ImplicitNewmarkSecondOrderTimeIntegrationRule time_rule_;
};
 
}  // namespace smith