state_variable_system.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/differentiable_numerics/field_store.hpp"
#include "smith/differentiable_numerics/nonlinear_block_solver.hpp"
#include "smith/differentiable_numerics/dirichlet_boundary_conditions.hpp"
#include "smith/differentiable_numerics/state_advancer.hpp"
#include "smith/differentiable_numerics/multiphysics_time_integrator.hpp"
#include "smith/differentiable_numerics/time_integration_rule.hpp"
#include "smith/physics/functional_weak_form.hpp"
#include "smith/differentiable_numerics/differentiable_physics.hpp"
#include "smith/physics/weak_form.hpp"
#include "smith/differentiable_numerics/system_base.hpp"
#include "smith/differentiable_numerics/coupling_params.hpp"
 
namespace smith {
 
template <int dim, typename StateSpace, typename InternalVarTimeRule = BackwardEulerFirstOrderTimeIntegrationRule,
          typename Coupling = std::tuple<>>
struct InternalVariableSystem : public SystemBase {
  using SystemBase::SystemBase;
 
  static_assert(InternalVarTimeRule::num_states == 2,
                "InternalVariableSystem requires a 2-state time integration rule");
 
  using InternalVariableWeakFormType =
      FunctionalWeakForm<dim, StateSpace, detail::TimeRuleParams<InternalVarTimeRule, StateSpace, Coupling>>;
 
  std::shared_ptr<InternalVariableWeakFormType> internal_variable_weak_form;  
  std::shared_ptr<DirichletBoundaryConditions> internal_variable_bc;          
  std::shared_ptr<InternalVarTimeRule> internal_variable_time_rule;           
  std::shared_ptr<const Coupling> coupling;                                   
 
  template <typename EvolutionType>
  void addEvolution(const std::string& domain_name, EvolutionType evolution_law)
  {
    auto captured_time_rule = internal_variable_time_rule;
    auto captured_coupling = coupling;
    internal_variable_weak_form->addBodyIntegral(domain_name, [=](auto t_info, auto /*X*/, auto... raw_args) {
      return detail::applyTimeRuleAndCoupling(
          *captured_time_rule, *captured_coupling, t_info,
          [&](auto alpha_current, auto alpha_dot, auto... interpolated_params) {
            auto residual_val =
                evolution_law(t_info, get<VALUE>(alpha_current), get<VALUE>(alpha_dot), interpolated_params...);
            tensor<double, dim> flux{};
            return smith::tuple{residual_val, flux};
          },
          raw_args...);
    });
  }
};
 
// ---------------------------------------------------------------------------
// Phase 1: registerInternalVariableFields
// ---------------------------------------------------------------------------
 
template <int dim, typename StateSpace, typename InternalVarTimeRule, typename... parameter_space>
auto registerInternalVariableFields(std::shared_ptr<FieldStore> field_store,
                                    FieldType<parameter_space>... parameter_types)
{
  auto internal_variable_time_rule = std::make_shared<InternalVarTimeRule>();
  FieldType<StateSpace> state_type("state_solve_state");
  field_store->addIndependent(state_type, internal_variable_time_rule);
  field_store->addDependent(state_type, FieldStore::TimeDerivative::VAL, "state");
 
  if constexpr (sizeof...(parameter_space) > 0) {
    auto prefix_param = [&](auto& pt) {
      pt.name = "param_" + pt.name;
      field_store->addParameter(pt);
    };
    (prefix_param(parameter_types), ...);
  }
 
  return PhysicsFields<dim, StateSpace::order, InternalVarTimeRule, StateSpace, StateSpace>{
      field_store, FieldType<StateSpace>(field_store->prefix("state_solve_state")),
      FieldType<StateSpace>(field_store->prefix("state"))};
}
 
template <int dim, typename StateSpace, typename InternalVarTimeRule, typename... parameter_space>
auto registerStateVariableFields(std::shared_ptr<FieldStore> field_store, FieldType<parameter_space>... parameter_types)
{
  return registerInternalVariableFields<dim, StateSpace, InternalVarTimeRule>(field_store, parameter_types...);
}
 
// ---------------------------------------------------------------------------
// Phase 2: buildInternalVariableSystem
// ---------------------------------------------------------------------------
 
namespace detail {
 
template <int dim, typename StateSpace, typename InternalVarTimeRule, typename Coupling>
  requires detail::is_coupling_packs_v<Coupling>
auto buildInternalVariableSystemImpl(std::shared_ptr<FieldStore> field_store, const Coupling& coupling,
                                     std::shared_ptr<SystemSolver> solver)
{
  auto internal_variable_time_rule = std::make_shared<InternalVarTimeRule>();
 
  FieldType<StateSpace> state_type(field_store->prefix("state_solve_state"), true);
  FieldType<StateSpace> state_old_type(field_store->prefix("state"));
 
  auto internal_variable_bc = field_store->getBoundaryConditions(state_type.name);
 
  using SystemType = InternalVariableSystem<dim, StateSpace, InternalVarTimeRule, Coupling>;
 
  std::string internal_variable_residual_name = field_store->prefix("state_residual");
  auto internal_variable_weak_form =
      detail::buildWeakFormWithCoupling<typename SystemType::InternalVariableWeakFormType>(
          field_store, internal_variable_residual_name, state_type.name, state_type, state_old_type,
          detail::flattenCouplingFields(coupling));
 
  auto sys = std::make_shared<SystemType>(field_store, solver,
                                          std::vector<std::shared_ptr<WeakForm>>{internal_variable_weak_form});
  sys->internal_variable_bc = internal_variable_bc;
  sys->internal_variable_time_rule = internal_variable_time_rule;
  sys->coupling = std::make_shared<Coupling>(coupling);
  sys->internal_variable_weak_form = internal_variable_weak_form;
 
  return sys;
}
 
}  // namespace detail
 
template <typename SelfFields>
  requires(detail::has_time_rule_v<SelfFields>)
auto buildInternalVariableSystem(std::shared_ptr<SystemSolver> solver, const SelfFields& self_fields)
{
  constexpr int dim = SelfFields::dim;
  using StateSpace = typename std::tuple_element_t<0, decltype(self_fields.fields)>::space_type;
  using InternalVarTimeRule = typename std::decay_t<SelfFields>::time_rule_type;
  auto field_store = self_fields.field_store;
  auto coupling = detail::collectCouplingFields();
  return detail::buildInternalVariableSystemImpl<dim, StateSpace, InternalVarTimeRule>(field_store, coupling, solver);
}
 
template <typename SelfFields, typename... PFs>
  requires(detail::has_time_rule_v<SelfFields>)
auto buildInternalVariableSystem(std::shared_ptr<SystemSolver> solver, const SelfFields& self_fields,
                                 const CouplingFields<PFs...>& coupled)
{
  constexpr int dim = SelfFields::dim;
  using StateSpace = typename std::tuple_element_t<0, decltype(self_fields.fields)>::space_type;
  using InternalVarTimeRule = typename std::decay_t<SelfFields>::time_rule_type;
  auto field_store = self_fields.field_store;
  auto coupling = detail::collectCouplingFields(coupled);
  return detail::buildInternalVariableSystemImpl<dim, StateSpace, InternalVarTimeRule>(field_store, coupling, solver);
}
 
template <typename SelfFields, typename... ParamSpaces>
  requires(detail::has_time_rule_v<SelfFields>)
auto buildInternalVariableSystem(std::shared_ptr<SystemSolver> solver, const SelfFields& self_fields,
                                 const ParamFields<ParamSpaces...>& params)
{
  constexpr int dim = SelfFields::dim;
  using StateSpace = typename std::tuple_element_t<0, decltype(self_fields.fields)>::space_type;
  using InternalVarTimeRule = typename std::decay_t<SelfFields>::time_rule_type;
  auto field_store = self_fields.field_store;
  auto coupling = detail::collectCouplingFields(params);
  return detail::buildInternalVariableSystemImpl<dim, StateSpace, InternalVarTimeRule>(field_store, coupling, solver);
}
 
template <typename SelfFields, typename... PFs, typename... ParamSpaces>
  requires(detail::has_time_rule_v<SelfFields>)
auto buildInternalVariableSystem(std::shared_ptr<SystemSolver> solver, const SelfFields& self_fields,
                                 const CouplingFields<PFs...>& coupled, const ParamFields<ParamSpaces...>& params)
{
  constexpr int dim = SelfFields::dim;
  using StateSpace = typename std::tuple_element_t<0, decltype(self_fields.fields)>::space_type;
  using InternalVarTimeRule = typename std::decay_t<SelfFields>::time_rule_type;
  auto field_store = self_fields.field_store;
  auto coupling = detail::collectCouplingFields(coupled, params);
  return detail::buildInternalVariableSystemImpl<dim, StateSpace, InternalVarTimeRule>(field_store, coupling, solver);
}
 
}  // namespace smith