thermomechanics.hpp

Go to the documentation of this file.

// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
// other Serac Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include "mfem.hpp"
 
#include "serac/physics/base_physics.hpp"
#include "serac/physics/thermomechanics_input.hpp"
#include "serac/physics/solid_mechanics.hpp"
#include "serac/physics/heat_transfer.hpp"
#include "serac/physics/materials/thermal_material.hpp"
#include "serac/physics/materials/solid_material.hpp"
 
namespace serac {
 
template <int order, int dim, typename... parameter_space>
class Thermomechanics : public BasePhysics {
public:
  Thermomechanics(const NonlinearSolverOptions thermal_nonlin_opts, const LinearSolverOptions thermal_lin_opts,
                  TimesteppingOptions thermal_timestepping, const NonlinearSolverOptions solid_nonlin_opts,
                  const LinearSolverOptions solid_lin_opts, TimesteppingOptions solid_timestepping,
                  GeometricNonlinearities geom_nonlin, const std::string& physics_name, std::string mesh_tag,
                  int cycle = 0, double time = 0.0)
      : Thermomechanics(
            std::make_unique<EquationSolver>(thermal_nonlin_opts, thermal_lin_opts,
                                             StateManager::mesh(mesh_tag).GetComm()),
            thermal_timestepping,
            std::make_unique<EquationSolver>(solid_nonlin_opts, solid_lin_opts, StateManager::mesh(mesh_tag).GetComm()),
            solid_timestepping, geom_nonlin, physics_name, mesh_tag, cycle, time)
  {
  }
 
  Thermomechanics(std::unique_ptr<EquationSolver> thermal_solver, TimesteppingOptions thermal_timestepping,
                  std::unique_ptr<EquationSolver> solid_solver, TimesteppingOptions solid_timestepping,
                  GeometricNonlinearities geom_nonlin, const std::string& physics_name, std::string mesh_tag,
                  int cycle = 0, double time = 0.0)
      : BasePhysics(physics_name, mesh_tag),
        thermal_(std::move(thermal_solver), thermal_timestepping, physics_name + "thermal", mesh_tag, {"displacement"},
                 cycle, time),
        solid_(std::move(solid_solver), solid_timestepping, geom_nonlin, physics_name + "mechanical", mesh_tag,
               {"temperature"}, cycle, time)
  {
    SLIC_ERROR_ROOT_IF(mesh_.Dimension() != dim,
                       axom::fmt::format("Compile time dimension and runtime mesh dimension mismatch"));
 
    states_.push_back(&thermal_.temperature());
    states_.push_back(&solid_.velocity());
    states_.push_back(&solid_.displacement());
  }
 
  Thermomechanics(const HeatTransferInputOptions& thermal_options, const SolidMechanicsInputOptions& solid_options,
                  const std::string& physics_name, std::string mesh_tag, int cycle = 0, double time = 0.0)
      : Thermomechanics(thermal_options.nonlin_solver_options, thermal_options.lin_solver_options,
                        thermal_options.timestepping_options, solid_options.nonlin_solver_options,
                        solid_options.lin_solver_options, solid_options.timestepping_options, solid_options.geom_nonlin,
                        physics_name, mesh_tag, cycle, time)
  {
  }
 
  Thermomechanics(const ThermomechanicsInputOptions& options, const std::string& physics_name, std::string mesh_tag,
                  int cycle = 0, double time = 0.0)
      : Thermomechanics(options.thermal_options, options.solid_options, physics_name, mesh_tag, cycle, time)
  {
    if (options.coef_thermal_expansion) {
      std::unique_ptr<mfem::Coefficient> cte(options.coef_thermal_expansion->constructScalar());
      std::unique_ptr<mfem::Coefficient> ref_temp(options.reference_temperature->constructScalar());
 
      // setThermalExpansion(std::move(cte), std::move(ref_temp));
    }
  }
 
  void completeSetup() override
  {
    thermal_.completeSetup();
    solid_.completeSetup();
  }
 
  void resetStates(int cycle = 0, double time = 0.0) override
  {
    BasePhysics::initializeBasePhysicsStates(cycle, time);
    thermal_.resetStates(cycle, time);
    solid_.resetStates(cycle, time);
  }
 
  const FiniteElementState& state(const std::string& state_name) const override
  {
    if (state_name == "displacement") {
      return solid_.displacement();
    } else if (state_name == "velocity") {
      return solid_.velocity();
    } else if (state_name == "temperature") {
      return thermal_.temperature();
    }
 
    SLIC_ERROR_ROOT(axom::fmt::format("State '{}' requested from solid mechanics module '{}', but it doesn't exist",
                                      state_name, name_));
    return solid_.displacement();
  }
 
  void setState(const std::string& state_name, const FiniteElementState& state) override
  {
    if (state_name == "displacement") {
      const_cast<FiniteElementState&>(solid_.displacement()) = state;
      return;
    } else if (state_name == "velocity") {
      const_cast<FiniteElementState&>(solid_.velocity()) = state;
      return;
    } else if (state_name == "temperature") {
      const_cast<FiniteElementState&>(thermal_.temperature()) = state;
      return;
    }
 
    SLIC_ERROR_ROOT(axom::fmt::format(
        "setState for state named '{}' requested from thermomechanics module '{}', but it doesn't exist", state_name,
        name_));
  }
 
  virtual std::vector<std::string> stateNames() const override
  {
    return std::vector<std::string>{"displacement", "velocity", "temperature"};
  }
 
  const FiniteElementState& adjoint(const std::string& state_name) const override
  {
    if (state_name == "displacement") {
      return solid_.adjoint("displacement");
    } else if (state_name == "temperature") {
      return thermal_.adjoint("temperature");
    }
 
    SLIC_ERROR_ROOT(axom::fmt::format("Adjoint '{}' requested from solid mechanics module '{}', but it doesn't exist",
                                      state_name, name_));
    return solid_.displacement();
  }
 
  void advanceTimestep(double dt) override
  {
    thermal_.setParameter(0, solid_.displacement());
    thermal_.advanceTimestep(dt);
 
    solid_.setParameter(0, thermal_.temperature());
    solid_.advanceTimestep(dt);
 
    cycle_ += 1;
    time_ += dt;
  }
 
  template <typename T>
  std::shared_ptr<QuadratureData<T>> createQuadratureDataBuffer(T initial_state)
  {
    return solid_.createQuadratureDataBuffer(initial_state);
  }
 
  template <typename ThermalMechanicalMaterial>
  struct ThermalMaterialInterface {
    using State = typename ThermalMechanicalMaterial::State;  
 
    const ThermalMechanicalMaterial mat;  
 
    ThermalMaterialInterface(const ThermalMechanicalMaterial& m) : mat(m)
    {
      // empty
    }
 
    template <typename T1, typename T2, typename T3, typename T4, typename... param_types>
    SERAC_HOST_DEVICE auto operator()(const T1& /* x */, const T2& temperature, const T3& temperature_gradient,
                                      const T4& displacement, param_types... parameters) const
    {
      // BT: this will not update the state correctly. I just want to get the code compiling before plumbing the
      // state variables.
      State state{};
 
      auto [u, du_dX]                 = displacement;
      auto [T, heat_capacity, s0, q0] = mat(state, du_dX, temperature, temperature_gradient, parameters...);
 
      return serac::tuple{heat_capacity, q0};
    }
  };
 
  template <typename ThermalMechanicalMaterial>
  struct MechanicalMaterialInterface {
    using State = typename ThermalMechanicalMaterial::State;  
 
    const ThermalMechanicalMaterial mat;  
 
    const double density;  
 
    MechanicalMaterialInterface(const ThermalMechanicalMaterial& m) : mat(m), density(m.density)
    {
      // empty
    }
 
    template <typename T1, typename T2, typename... param_types>
    SERAC_HOST_DEVICE auto operator()(State& state, const T1& displacement_gradient, const T2& temperature,
                                      param_types... parameters) const
    {
      auto [theta, dtheta_dX]         = temperature;
      auto [T, heat_capacity, s0, q0] = mat(state, displacement_gradient, theta, dtheta_dX, parameters...);
      return T;
    }
  };
 
  template <int... active_parameters, typename MaterialType, typename StateType>
  void setMaterial(DependsOn<active_parameters...>, MaterialType material,
                   std::shared_ptr<QuadratureData<StateType>> qdata)
  {
    // note: these parameter indices are offset by 1 since, internally, this module uses the first parameter
    // to communicate the temperature and displacement field information to the other physics module
    //
    thermal_.setMaterial(DependsOn<0, active_parameters + 1 ...>{}, ThermalMaterialInterface<MaterialType>{material});
    solid_.setMaterial(DependsOn<0, active_parameters + 1 ...>{}, MechanicalMaterialInterface<MaterialType>{material},
                       qdata);
  }
 
  template <typename MaterialType, typename StateType = Empty>
  void setMaterial(MaterialType material, std::shared_ptr<QuadratureData<StateType>> qdata = EmptyQData)
  {
    setMaterial(DependsOn<>{}, material, qdata);
  }
 
  void setTemperatureBCs(const std::set<int>&                                   temperature_attributes,
                         std::function<double(const mfem::Vector& x, double t)> prescribed_value)
  {
    thermal_.setTemperatureBCs(temperature_attributes, prescribed_value);
  }
 
  void setDisplacementBCs(const std::set<int>&                                           displacement_attributes,
                          std::function<void(const mfem::Vector& x, mfem::Vector& disp)> prescribed_value)
  {
    solid_.setDisplacementBCs(displacement_attributes, prescribed_value);
  }
 
  template <typename FluxType>
  void setHeatFluxBCs(FluxType flux_function)
  {
    thermal_.setFluxBCs(flux_function);
  }
 
  void setDisplacement(std::function<void(const mfem::Vector& x, mfem::Vector& u)> displacement)
  {
    solid_.setDisplacement(displacement);
  }
 
  void setTemperature(std::function<double(const mfem::Vector& x, double t)> temperature)
  {
    thermal_.setTemperature(temperature);
  }
 
  template <typename BodyForceType>
  void addBodyForce(BodyForceType body_force_function)
  {
    solid_.addBodyForce(body_force_function);
  }
 
  template <typename HeatSourceType>
  void addHeatSource(HeatSourceType source_function)
  {
    thermal_.setSource(source_function);
  }
 
  const serac::FiniteElementState& displacement() const { return solid_.displacement(); };
 
  const serac::FiniteElementState& temperature() const { return thermal_.temperature(); };
 
protected:
  using displacement_field = H1<order, dim>;  
  using temperature_field  = H1<order>;       
 
  HeatTransfer<order, dim, Parameters<displacement_field, parameter_space...>> thermal_;
 
  SolidMechanics<order, dim, Parameters<temperature_field, parameter_space...>> solid_;
};
 
}  // namespace serac