thermal_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/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, int temp_order, typename TemperatureTimeRule = QuasiStaticFirstOrderTimeIntegrationRule,
          typename Coupling = std::tuple<>>
struct ThermalSystem : public SystemBase {
  using SystemBase::SystemBase;
 
  static_assert(TemperatureTimeRule::num_states == 2, "ThermalSystem requires a 2-state time integration rule");
 
  using ThermalWeakFormType =
      FunctionalWeakForm<dim, H1<temp_order>, detail::TimeRuleParams<TemperatureTimeRule, H1<temp_order>, Coupling>>;
 
  std::shared_ptr<ThermalWeakFormType> thermal_weak_form;       
  std::shared_ptr<DirichletBoundaryConditions> temperature_bc;  
  std::shared_ptr<TemperatureTimeRule> temperature_time_rule;   
  std::shared_ptr<const Coupling> coupling;                     
 
  template <typename MaterialType>
  void setMaterial(const MaterialType& material, const std::string& domain_name)
  {
    auto captured_temp_rule = temperature_time_rule;
    auto captured_coupling = coupling;
 
    thermal_weak_form->addBodyIntegral(domain_name, [=](auto t_info, auto /*X*/, auto... raw_args) {
      return detail::applyTimeRuleAndCoupling(
          *captured_temp_rule, *captured_coupling, t_info,
          [&](auto T_current, auto T_dot, auto... interpolated_params) {
            auto [heat_capacity, heat_flux] =
                material(t_info, get<VALUE>(T_current), get<DERIVATIVE>(T_current), interpolated_params...);
            return smith::tuple{heat_capacity * get<VALUE>(T_dot), -heat_flux};
          },
          raw_args...);
    });
  }
 
  template <typename MaterialType>
  void setMaterialAndHeatSource(const MaterialType& material, const std::string& domain_name)
  {
    auto captured_temp_rule = temperature_time_rule;
    auto captured_coupling = coupling;
 
    thermal_weak_form->addBodyIntegral(domain_name, [=](auto t_info, auto /*X*/, auto... raw_args) {
      return detail::applyTimeRuleAndCoupling(
          *captured_temp_rule, *captured_coupling, t_info,
          [&](auto T_current, auto T_dot, auto... interpolated_params) {
            auto [heat_capacity, heat_flux, heat_source] =
                material(t_info, get<VALUE>(T_current), get<DERIVATIVE>(T_current), interpolated_params...);
            return smith::tuple{heat_capacity * get<VALUE>(T_dot) - heat_source, -heat_flux};
          },
          raw_args...);
    });
  }
 
  template <typename HeatSourceType>
  void addHeatSource(const std::string& domain_name, HeatSourceType source_function)
  {
    auto captured_rule = temperature_time_rule;
    auto captured_coupling = coupling;
    thermal_weak_form->addBodySource(domain_name, [=](auto t_info, auto X, auto... raw_args) {
      return detail::applyTimeRuleAndCoupling(
          *captured_rule, *captured_coupling, t_info,
          [&](auto T_current, auto /*T_dot*/, auto... interpolated_params) {
            return source_function(t_info.time(), X, T_current, interpolated_params...);
          },
          raw_args...);
    });
  }
 
  template <typename HeatFluxType>
  void addHeatFlux(const std::string& boundary_name, HeatFluxType flux_function)
  {
    auto captured_rule = temperature_time_rule;
    auto captured_coupling = coupling;
    thermal_weak_form->addBoundaryFlux(boundary_name, [=](auto t_info, auto X, auto n, auto... raw_args) {
      return detail::applyTimeRuleAndCoupling(
          *captured_rule, *captured_coupling, t_info,
          [&](auto T_current, auto /*T_dot*/, auto... interpolated_params) {
            return flux_function(t_info.time(), X, n, T_current, interpolated_params...);
          },
          raw_args...);
    });
  }
 
  void setTemperatureBC(const Domain& domain) { temperature_bc->template setFixedScalarBCs<dim>(domain); }
 
  template <typename AppliedTemperatureFunction>
  void setTemperatureBC(const Domain& domain, AppliedTemperatureFunction f)
  {
    temperature_bc->template setScalarBCs<dim>(domain, f);
  }
 
 private:
};
 
struct ThermalOptions {};
 
template <int dim, int temp_order, typename TemperatureTimeRule>
auto registerThermalFields(std::shared_ptr<FieldStore> field_store,
                           const ThermalOptions& /*options*/ = ThermalOptions{})
{
  FieldType<H1<1, dim>> shape_disp_type("shape_displacement");
  if (!field_store->hasField(field_store->prefix(shape_disp_type.name))) {
    field_store->addShapeDisp(shape_disp_type);
  }
 
  auto temperature_time_rule_ptr = std::make_shared<TemperatureTimeRule>();
  FieldType<H1<temp_order>> temperature_type("temperature_solve_state");
  field_store->addIndependent(temperature_type, temperature_time_rule_ptr);
  field_store->addDependent(temperature_type, FieldStore::TimeDerivative::VAL, "temperature");
 
  return PhysicsFields<dim, temp_order, TemperatureTimeRule, H1<temp_order>, H1<temp_order>>{
      field_store, FieldType<H1<temp_order>>(field_store->prefix("temperature_solve_state")),
      FieldType<H1<temp_order>>(field_store->prefix("temperature"))};
}
 
namespace detail {
 
template <int dim, int temp_order, typename TemperatureTimeRule, typename Coupling>
  requires detail::is_coupling_packs_v<Coupling>
auto buildThermalSystemImpl(std::shared_ptr<FieldStore> field_store, const Coupling& coupling,
                            std::shared_ptr<SystemSolver> solver, const ThermalOptions& /*options*/)
{
  auto temperature_time_rule_ptr = std::make_shared<TemperatureTimeRule>();
 
  FieldType<H1<1, dim>> shape_disp_type(field_store->prefix("shape_displacement"));
  FieldType<H1<temp_order>> temperature_type(field_store->prefix("temperature_solve_state"), true);
  FieldType<H1<temp_order>> temperature_old_type(field_store->prefix("temperature"));
 
  auto temperature_bc = field_store->getBoundaryConditions(temperature_type.name);
 
  using SystemType = ThermalSystem<dim, temp_order, TemperatureTimeRule, Coupling>;
 
  std::string thermal_flux_name = field_store->prefix("thermal_flux");
  auto thermal_weak_form = detail::buildWeakFormWithCoupling<typename SystemType::ThermalWeakFormType>(
      field_store, thermal_flux_name, temperature_type.name, temperature_type, temperature_old_type,
      detail::flattenCouplingFields(coupling));
 
  auto sys =
      std::make_shared<SystemType>(field_store, solver, std::vector<std::shared_ptr<WeakForm>>{thermal_weak_form});
  sys->temperature_bc = temperature_bc;
  sys->temperature_time_rule = temperature_time_rule_ptr;
  sys->coupling = std::make_shared<Coupling>(coupling);
  sys->thermal_weak_form = thermal_weak_form;
 
  return sys;
}
 
}  // namespace detail
 
template <typename SelfFields>
  requires(detail::has_time_rule_v<SelfFields>)
auto buildThermalSystem(std::shared_ptr<SystemSolver> solver, const ThermalOptions& options,
                        const SelfFields& self_fields)
{
  constexpr int dim = SelfFields::dim;
  constexpr int temp_order = SelfFields::order;
  using TemperatureTimeRule = typename std::decay_t<SelfFields>::time_rule_type;
  auto field_store = self_fields.field_store;
  auto coupling = detail::collectCouplingFields();
  return detail::buildThermalSystemImpl<dim, temp_order, TemperatureTimeRule>(field_store, coupling, solver, options);
}
 
template <typename SelfFields, typename... PFs>
  requires(detail::has_time_rule_v<SelfFields>)
auto buildThermalSystem(std::shared_ptr<SystemSolver> solver, const ThermalOptions& options,
                        const SelfFields& self_fields, const CouplingFields<PFs...>& coupled)
{
  constexpr int dim = SelfFields::dim;
  constexpr int temp_order = SelfFields::order;
  using TemperatureTimeRule = typename std::decay_t<SelfFields>::time_rule_type;
  auto field_store = self_fields.field_store;
  auto coupling = detail::collectCouplingFields(coupled);
  return detail::buildThermalSystemImpl<dim, temp_order, TemperatureTimeRule>(field_store, coupling, solver, options);
}
 
template <typename SelfFields, typename... ParamSpaces>
  requires(detail::has_time_rule_v<SelfFields>)
auto buildThermalSystem(std::shared_ptr<SystemSolver> solver, const ThermalOptions& options,
                        const SelfFields& self_fields, const ParamFields<ParamSpaces...>& params)
{
  constexpr int dim = SelfFields::dim;
  constexpr int temp_order = SelfFields::order;
  using TemperatureTimeRule = typename std::decay_t<SelfFields>::time_rule_type;
  auto field_store = self_fields.field_store;
  auto coupling = detail::collectCouplingFields(params);
  return detail::buildThermalSystemImpl<dim, temp_order, TemperatureTimeRule>(field_store, coupling, solver, options);
}
 
template <typename SelfFields, typename... PFs, typename... ParamSpaces>
  requires(detail::has_time_rule_v<SelfFields>)
auto buildThermalSystem(std::shared_ptr<SystemSolver> solver, const ThermalOptions& options,
                        const SelfFields& self_fields, const CouplingFields<PFs...>& coupled,
                        const ParamFields<ParamSpaces...>& params)
{
  constexpr int dim = SelfFields::dim;
  constexpr int temp_order = SelfFields::order;
  using TemperatureTimeRule = typename std::decay_t<SelfFields>::time_rule_type;
  auto field_store = self_fields.field_store;
  auto coupling = detail::collectCouplingFields(coupled, params);
  return detail::buildThermalSystemImpl<dim, temp_order, TemperatureTimeRule>(field_store, coupling, solver, options);
}
 
template <int dim, int temp_order, typename TemperatureTimeRule = QuasiStaticFirstOrderTimeIntegrationRule,
          typename... ParamSpaces>
auto buildThermalSystem(const NonlinearSolverOptions& nonlinear_opts, const LinearSolverOptions& linear_opts,
                        const ThermalOptions& options, std::shared_ptr<smith::Mesh> mesh,
                        FieldType<ParamSpaces>... params)
{
  auto field_store = std::make_shared<FieldStore>(mesh);
  auto solver = std::make_shared<SystemSolver>(buildNonlinearBlockSolver(nonlinear_opts, linear_opts, *mesh));
  auto thermal_fields = registerThermalFields<dim, temp_order, TemperatureTimeRule>(field_store, options);
  if constexpr (sizeof...(ParamSpaces) > 0) {
    auto param_fields = registerParameterFields(field_store, std::move(params)...);
    return buildThermalSystem(solver, options, thermal_fields, param_fields);
  } else {
    return buildThermalSystem(solver, options, thermal_fields);
  }
}
 
}  // namespace smith