equation_of_state.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/vectorization.h>
 
#include <meltpooldg/compressible_flow/material.hpp>
#include <meltpooldg/utilities/concepts.hpp>
#include <meltpooldg/utilities/dealii_tensor.hpp>
 
#include <cmath>
#include <utility>
 
namespace MeltPoolDG::Flow
{
  template <typename T>
  concept EOSIsValueView = requires(const T v) {
    {
      v.density()
    };
    {
      v.velocity()
    };
    {
      v.total_energy()
    };
  };
 
  template <typename T>
  concept EOSIsGradientView = requires(const T g) {
    {
      g.grad_density()
    };
    {
      g.grad_velocity()
    };
    {
      g.grad_total_energy()
    };
  };
 
  template <typename T>
  concept IdealGasIsMaterialView = requires(const T m) {
    {
      m.heat_capacity_ratio()
    };
    {
      m.specific_gas_constant()
    };
  };
 
  template <typename T>
  concept StiffenedGasIsMaterialView = requires(const T m) {
    {
      m.heat_capacity_ratio()
    };
    {
      m.specific_isobaric_heat()
    };
    {
      m.stiffening_pressure()
    };
  };
 
  template <typename T>
  concept NobleAbelStiffenedGasIsMaterialView = requires(const T m) {
    {
      m.heat_capacity_ratio()
    };
    {
      m.specific_isobaric_heat()
    };
    {
      m.stiffening_pressure()
    };
    {
      m.heat_bound()
    };
    {
      m.covolume()
    };
  };
 
  struct IdealGasEOS
  {
    template <EOSIsValueView ValueView, IdealGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      thermodynamic_pressure(const ValueView &value_view, const MaterialView &material_view)
    {
      return (material_view.heat_capacity_ratio() - 1.) *
             (value_view.total_energy() -
              value_view.density() * 0.5 *
                scalar_product(value_view.velocity(), value_view.velocity()));
    }
 
    template <EOSIsValueView         ValueView,
              EOSIsGradientView      GradientView,
              IdealGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      grad_temperature(const ValueView    &value_view,
                       const GradientView &gradient_view,
                       const MaterialView &material_view)
    {
      const auto grad_E =
        1. / value_view.density() *
        (gradient_view.grad_total_energy() -
         1. / value_view.density() * value_view.total_energy() * gradient_view.grad_density());
 
      return (material_view.heat_capacity_ratio() - 1.0) / material_view.specific_gas_constant() *
             (grad_E - matrix_vector_product(gradient_view.grad_velocity(), value_view.velocity()));
    }
 
    template <EOSIsValueView ValueView, IdealGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      speed_of_sound(const ValueView &value_view, const MaterialView &material_view)
    {
      return std::sqrt(material_view.heat_capacity_ratio() *
                       thermodynamic_pressure(value_view, material_view) / value_view.density());
    }
 
    template <EOSIsValueView ValueView, IdealGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      temperature(const ValueView &value_view, const MaterialView &material_view)
    {
      return thermodynamic_pressure(value_view, material_view) /
             (material_view.specific_gas_constant() * value_view.density());
    }
 
    template <typename ValueType, EOSIsValueView ValueView, IdealGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      inner_energy_from_pressure(const ValueType &pressure,
                                 const ValueView &,
                                 const MaterialView &material_view)
    {
      return pressure / (material_view.heat_capacity_ratio() - 1.);
    }
 
    template <typename ValueView, IdealGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      specific_inner_energy(const ValueView &value_view, const MaterialView &material_view)
    {
      return thermodynamic_pressure(value_view, material_view) /
             (value_view.density() * (material_view.heat_capacity_ratio() - 1.));
    }
  };
 
  struct StiffenedGasEOS
  {
    template <EOSIsValueView ValueView, StiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      thermodynamic_pressure(const ValueView &value_view, const MaterialView &material_view)
    {
      return (material_view.heat_capacity_ratio() - 1.) *
               (value_view.total_energy() -
                value_view.density() * 0.5 *
                  scalar_product(value_view.velocity(), value_view.velocity())) -
             material_view.heat_capacity_ratio() * material_view.stiffening_pressure();
    }
 
    template <EOSIsValueView             ValueView,
              EOSIsGradientView          GradientView,
              StiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      grad_temperature(const ValueView    &value_view,
                       const GradientView &gradient_view,
                       const MaterialView &material_view)
    {
      const auto inv_rho = 1. / value_view.density();
 
      const auto grad_E =
        inv_rho * (gradient_view.grad_total_energy() -
                   inv_rho * value_view.total_energy() * gradient_view.grad_density());
 
      return material_view.heat_capacity_ratio() / material_view.specific_isobaric_heat() *
             (grad_E - matrix_vector_product(gradient_view.grad_velocity(), value_view.velocity()) +
              material_view.stiffening_pressure() * inv_rho * inv_rho *
                gradient_view.grad_density());
    }
 
    template <EOSIsValueView ValueView, StiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      speed_of_sound(const ValueView &value_view, const MaterialView &material_view)
    {
      return std::sqrt(
        material_view.heat_capacity_ratio() *
        (thermodynamic_pressure(value_view, material_view) + material_view.stiffening_pressure()) /
        value_view.density());
    }
 
    template <EOSIsValueView ValueView, StiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      temperature(const ValueView &value_view, const MaterialView &material_view)
    {
      return (thermodynamic_pressure(value_view, material_view) +
              material_view.stiffening_pressure()) *
             material_view.heat_capacity_ratio() /
             (material_view.specific_isobaric_heat() * value_view.density() *
              (material_view.heat_capacity_ratio() - 1.));
    }
 
    template <typename ValueType, EOSIsValueView ValueView, StiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      inner_energy_from_pressure(const ValueType &pressure,
                                 const ValueView &,
                                 const MaterialView &material_view)
    {
      return (pressure +
              material_view.heat_capacity_ratio() * material_view.stiffening_pressure()) /
             (material_view.heat_capacity_ratio() - 1.);
    }
 
    template <typename ValueView, StiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      specific_inner_energy(const ValueView &value_view, const MaterialView &material_view)
    {
      return (thermodynamic_pressure(value_view, material_view) +
              material_view.heat_capacity_ratio() * material_view.stiffening_pressure()) /
             (value_view.density() * (material_view.heat_capacity_ratio() - 1.));
    }
  };
 
  struct NobleAbelStiffenedGasEOS
  {
    template <EOSIsValueView ValueView, NobleAbelStiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      thermodynamic_pressure(const ValueView &value_view, const MaterialView &material_view)
    {
      return ((material_view.heat_capacity_ratio() - 1.) *
                (value_view.total_energy() -
                 0.5 * value_view.density() *
                   scalar_product(value_view.velocity(), value_view.velocity()) -
                 value_view.density() * material_view.heat_bound()) -
              material_view.heat_capacity_ratio() * material_view.stiffening_pressure() *
                (1. - value_view.density() * material_view.covolume())) /
             (1. - value_view.density() * material_view.covolume());
    }
 
    template <EOSIsValueView                      ValueView,
              EOSIsGradientView                   GradientView,
              NobleAbelStiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      grad_temperature(const ValueView    &value_view,
                       const GradientView &gradient_view,
                       const MaterialView &material_view)
    {
      const auto inv_rho = 1. / value_view.density();
 
      const auto grad_E =
        inv_rho * (gradient_view.grad_total_energy() -
                   inv_rho * value_view.total_energy() * gradient_view.grad_density());
 
      return material_view.heat_capacity_ratio() / material_view.specific_isobaric_heat() *
             (grad_E - matrix_vector_product(gradient_view.grad_velocity(), value_view.velocity()) +
              material_view.stiffening_pressure() * inv_rho * inv_rho *
                gradient_view.grad_density());
    }
 
    template <EOSIsValueView ValueView, NobleAbelStiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      speed_of_sound(const ValueView &value_view, const MaterialView &material_view)
    {
      return std::sqrt(
        material_view.heat_capacity_ratio() *
        (thermodynamic_pressure(value_view, material_view) + material_view.stiffening_pressure()) /
        (value_view.density() * (1. - value_view.density() * material_view.covolume())));
    }
 
    template <EOSIsValueView ValueView, NobleAbelStiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      temperature(const ValueView &value_view, const MaterialView &material_view)
    {
      return (thermodynamic_pressure(value_view, material_view) +
              material_view.stiffening_pressure()) *
             (1. / value_view.density() - material_view.covolume()) *
             material_view.heat_capacity_ratio() / (material_view.heat_capacity_ratio() - 1.) *
             material_view.specific_isobaric_heat();
    }
 
    template <typename ValueType,
              EOSIsValueView                      ValueView,
              NobleAbelStiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      inner_energy_from_pressure(const ValueType    &pressure,
                                 const ValueView    &value_view,
                                 const MaterialView &material_view)
    {
      return (pressure +
              material_view.heat_capacity_ratio() * material_view.stiffening_pressure()) /
               (material_view.heat_capacity_ratio() - 1.) *
               (1. - value_view.density() * material_view.covolume()) +
             value_view.density() * material_view.heat_bound();
    }
 
    template <typename ValueView, NobleAbelStiffenedGasIsMaterialView MaterialView>
    static inline DEAL_II_ALWAYS_INLINE //
      auto
      specific_inner_energy(const ValueView &value_view, const MaterialView &material_view)
    {
      return (thermodynamic_pressure(value_view, material_view) +
              material_view.heat_capacity_ratio() * material_view.stiffening_pressure()) /
               (material_view.heat_capacity_ratio() - 1.) *
               (1. / value_view.density() - material_view.covolume()) +
             material_view.heat_bound();
    }
  };
 
  template <typename T>
  concept HasSupportedEOS = requires {
    {
      T::supported_eos
    } -> std::convertible_to<const std::array<CompressibleFlow::EquationOfState,
                                              std::tuple_size_v<decltype(T::supported_eos)>> &>;
  };
 
  template <typename Derived, CompressibleFlow::EquationOfState EOS>
  constexpr bool
  supports_eos()
  {
    return std::ranges::find(Derived::supported_eos, EOS) != Derived::supported_eos.end();
  }
 
  template <typename Derived, typename F>
  inline DEAL_II_ALWAYS_INLINE decltype(auto)
  dispatch_eos(CompressibleFlow::EquationOfState eos_type, F &&f)
  {
    auto check_support = [](CompressibleFlow::EquationOfState type) {
      if constexpr (HasSupportedEOS<Derived>)
        {
          if (type == CompressibleFlow::EquationOfState::ideal_gas)
            return supports_eos<Derived, CompressibleFlow::EquationOfState::ideal_gas>();
          if (type == CompressibleFlow::EquationOfState::stiffened_gas)
            return supports_eos<Derived, CompressibleFlow::EquationOfState::stiffened_gas>();
          if (type == CompressibleFlow::EquationOfState::noble_abel_stiffened_gas)
            return supports_eos<Derived,
                                CompressibleFlow::EquationOfState::noble_abel_stiffened_gas>();
          return false;
        }
      return true;
    };
 
    switch (eos_type)
      {
        case CompressibleFlow::EquationOfState::ideal_gas:
          if constexpr (check_support(CompressibleFlow::EquationOfState::ideal_gas))
            return std::forward<F>(f)(IdealGasEOS{});
          break;
 
        case CompressibleFlow::EquationOfState::stiffened_gas:
          if constexpr (check_support(CompressibleFlow::EquationOfState::stiffened_gas))
            return std::forward<F>(f)(StiffenedGasEOS{});
          break;
 
        case CompressibleFlow::EquationOfState::noble_abel_stiffened_gas:
          if constexpr (check_support(CompressibleFlow::EquationOfState::noble_abel_stiffened_gas))
            return std::forward<F>(f)(NobleAbelStiffenedGasEOS{});
          break;
 
        default:
          break;
      }
 
    AssertThrow(false, dealii::ExcMessage("The provided EOS is not supported by this View."));
  }
 
 
 
  template <int dim, ArithmeticType ValueType, typename Derived>
  struct EOSValueMixin
  {
    decltype(auto)
    pressure() const
    {
      return dispatch_eos<Derived>(eos_type(), [&](auto eos) -> decltype(auto) {
        return eos.thermodynamic_pressure(derived(), derived());
      });
    }
 
    decltype(auto)
    temperature() const
    {
      return dispatch_eos<Derived>(eos_type(), [&](auto eos) -> decltype(auto) {
        return eos.temperature(derived(), derived());
      });
    }
 
    decltype(auto)
    speed_of_sound() const
    {
      return dispatch_eos<Derived>(eos_type(), [&](auto eos) -> decltype(auto) {
        return eos.speed_of_sound(derived(), derived());
      });
    }
 
    decltype(auto)
    inner_energy_from_pressure(const ValueType &pressure) const
    {
      return dispatch_eos<Derived>(eos_type(), [&](auto eos) -> decltype(auto) {
        return eos.inner_energy_from_pressure(pressure, derived(), derived());
      });
    }
 
    decltype(auto)
    specific_inner_energy() const
    {
      return dispatch_eos<Derived>(eos_type(), [&](auto eos) -> decltype(auto) {
        return eos.specific_inner_energy(derived(), derived());
      });
    }
 
  private:
    decltype(auto)
    eos_type() const
    {
      return static_cast<const Derived &>(*this).eos_type();
    }
 
    const Derived &
    derived() const
    {
      return static_cast<const Derived &>(*this);
    }
  };
 
  template <int dim, typename Derived>
  struct EOSGradientMixin
  {
    decltype(auto)
    grad_temperature() const
    {
      return dispatch_eos<Derived>(eos_type(), [&](auto eos) -> decltype(auto) {
        return eos.grad_temperature(derived(), derived(), derived());
      });
    }
 
  private:
    const Derived &
    derived() const
    {
      return static_cast<const Derived &>(*this);
    }
 
    decltype(auto)
    eos_type() const
    {
      return static_cast<const Derived &>(*this).eos_type();
    }
  };
} // namespace MeltPoolDG::Flow