delta_approximation_phase_weighted.hpp

Go to the documentation of this file.

#pragma once
#include <deal.II/base/vectorization.h>
 
#include <meltpooldg/level_set/delta_approximation_phase_weighted_data.hpp>
#include <meltpooldg/level_set/level_set_tools.hpp>
 
namespace MeltPoolDG::LevelSet
{
  template <typename number>
  class DeltaApproximationBase
  {
  public:
    virtual ~DeltaApproximationBase() = default;
 
    virtual number
    compute_weight(const number ls_heaviside) const = 0;
 
    virtual dealii::VectorizedArray<number>
    compute_weight(const dealii::VectorizedArray<number> &ls_heaviside) const = 0;
  };
 
 
  template <typename number>
  class DeltaApproximationHeavyPhaseOnly : public DeltaApproximationBase<number>
  {
  public:
    DeltaApproximationHeavyPhaseOnly() = default;
 
    inline number
    compute_weight(const number level_set_heaviside) const override
    {
      return compute_weight_internal(level_set_heaviside);
    }
 
    inline dealii::VectorizedArray<number>
    compute_weight(const dealii::VectorizedArray<number> &level_set_heaviside) const override
    {
      dealii::VectorizedArray<number> temp(0.0);
      for (unsigned int v = 0; v < dealii::VectorizedArray<number>::size(); ++v)
        temp = compute_weight_internal(level_set_heaviside[v]);
      return temp;
    }
 
  private:
    template <typename value_type>
    inline value_type
    compute_weight_internal(const value_type &level_set_heaviside) const
    {
      if (level_set_heaviside > 0.5)
        return 2.0;
      else
        return 0.0;
    }
  };
 
 
  template <typename number>
  class DeltaApproximationHeavisidePhaseWeighted : public DeltaApproximationBase<number>
  {
  public:
    DeltaApproximationHeavisidePhaseWeighted(
      const DeltaApproximationPhaseWeightedData<number> &data)
      : w_g(data.gas_phase_weight)
      , w_h(data.heavy_phase_weight)
      , correction_factor(2. / (w_g + w_h))
    {
      AssertThrow(w_g + w_h != 0.0,
                  dealii::ExcMessage("When using a phase weighted Dirac delta function "
                                     "approximation use weights that don't sum up to zero! "
                                     "Abort..."));
    }
 
    inline number
    compute_weight(const number level_set_heaviside) const override
    {
      return compute_weight_internal(level_set_heaviside);
    }
 
    inline dealii::VectorizedArray<number>
    compute_weight(const dealii::VectorizedArray<number> &level_set_heaviside) const override
    {
      return compute_weight_internal(level_set_heaviside);
    }
 
  private:
    template <typename value_type>
    inline value_type
    compute_weight_internal(const value_type &level_set_heaviside) const
    {
      return Tools::interpolate(level_set_heaviside, w_g, w_h) * correction_factor;
    }
 
    const number w_g;
    const number w_h;
 
    /*
     * correction factor
     *
     *      2
     * -----------
     *  w_g + w_h
     */
    const number correction_factor;
  };
 
  template <typename number>
  class DeltaApproximationHeavisideTimesHeavisidePhaseWeighted
    : public DeltaApproximationBase<number>
  {
  public:
    DeltaApproximationHeavisideTimesHeavisidePhaseWeighted(
      const DeltaApproximationPhaseWeightedData<number> &data)
      : w_1g(data.gas_phase_weight)
      , w_1h(data.heavy_phase_weight)
      , w_2g(data.gas_phase_weight_2)
      , w_2h(data.heavy_phase_weight_2)
      , correction_factor(6. / (2. * w_1g * w_2g + w_1g * w_2h + w_1h * w_2g + 2. * w_1h * w_2h))
    {
      AssertThrow(std::abs(2. * w_1g * w_2g + w_1g * w_2h + w_1h * w_2g + 2. * w_1h * w_2h) >
                    std::numeric_limits<number>::epsilon(),
                  dealii::ExcMessage("When using a phase quadratic weighted Dirac delta "
                                     "function approximation use weights that fulfill this "
                                     "condition! Abort..."));
    }
 
    inline number
    compute_weight(const number level_set_heaviside) const override
    {
      return compute_weight_internal(level_set_heaviside);
    }
 
    inline dealii::VectorizedArray<number>
    compute_weight(const dealii::VectorizedArray<number> &level_set_heaviside) const override
    {
      return compute_weight_internal(level_set_heaviside);
    }
 
  private:
    template <typename value_type>
    inline value_type
    compute_weight_internal(const value_type &level_set_heaviside) const
    {
      const auto w_1 = Tools::interpolate(level_set_heaviside, w_1g, w_1h);
      const auto w_2 = Tools::interpolate(level_set_heaviside, w_2g, w_2h);
      return w_1 * w_2 * correction_factor;
    }
 
    const number w_1g;
    const number w_1h;
    const number w_2g;
    const number w_2h;
 
    /*
     * correction factor
     *
     *                          6
     * ---------------------------------------------------
     *  2 w_1g w_2g + w_1g w_2h + w_1h w_2g + 2 w_1h w_2h
     */
    const number correction_factor;
  };
 
  template <typename number>
  class DeltaApproximationReciprocalPhaseWeighted : public DeltaApproximationBase<number>
  {
  public:
    DeltaApproximationReciprocalPhaseWeighted(
      const DeltaApproximationPhaseWeightedData<number> &data)
      : w_g(data.gas_phase_weight)
      , w_h(data.heavy_phase_weight)
      , correction_factor((w_g - w_h) / (w_g * w_h * std::log(w_g / w_h)))
    {
      AssertThrow(w_g > 0.0,
                  dealii::ExcMessage("When using the Dirac delta function approximation "
                                     "weighted consistent with evaporation use positive weights! "
                                     "Abort..."));
      AssertThrow(w_h > 0.0,
                  dealii::ExcMessage("When using the Dirac delta function approximation "
                                     "weighted consistent with evaporation use positive weights! "
                                     "Abort..."));
      AssertThrow(std::abs(w_g - w_h) > std::numeric_limits<number>::epsilon(),
                  dealii::ExcMessage("When using the Dirac delta function approximation "
                                     "weighted consistent with evaporation the phase weights must "
                                     "differ! Abort..."));
    }
 
    inline number
    compute_weight(const number level_set_heaviside) const override
    {
      return compute_weight_internal(level_set_heaviside);
    }
 
    inline dealii::VectorizedArray<number>
    compute_weight(const dealii::VectorizedArray<number> &level_set_heaviside) const override
    {
      return compute_weight_internal(level_set_heaviside);
    }
 
  private:
    template <typename value_type>
    inline value_type
    compute_weight_internal(const value_type &level_set_heaviside) const
    {
      return Tools::interpolate_reciprocal(level_set_heaviside, w_g, w_h) * correction_factor;
    }
 
    const number w_g;
    const number w_h;
 
    /*
     * correction factor
     *
     *         w_g - w_h
     * -------------------------
     *  w_g w_h ln( w_g / w_h )
     */
    const number correction_factor;
  };
 
  template <typename number>
  class DeltaApproximationReciprocalTimesHeavisidePhaseWeighted
    : public DeltaApproximationBase<number>
  {
  public:
    DeltaApproximationReciprocalTimesHeavisidePhaseWeighted(
      const DeltaApproximationPhaseWeightedData<number> &data)
      : w_1g(data.gas_phase_weight)
      , w_1h(data.heavy_phase_weight)
      , w_2g(data.gas_phase_weight_2)
      , w_2h(data.heavy_phase_weight_2)
      , correction_factor(
          // clang-format off
          1. / (
              w_2g * (w_1g * w_1h * std::log(w_1g / w_1h))
              / (w_1g - w_1h)
            +
              (w_2h - w_2g) * (1. / w_1g * (std::log(w_1h / w_1g) - 1.) + 1. / w_1h)
              / dealii::Utilities::fixed_power<2>(1. / w_1h - 1. / w_1g)
          )
          // clang-format on
        )
    {
      AssertThrow(w_1g > 0.0,
                  dealii::ExcMessage(
                    "When using the Dirac delta function approximation that is phase weighted"
                    "reciprocal times heaviside use positive weights! Abort..."));
      AssertThrow(w_1h > 0.0,
                  dealii::ExcMessage(
                    "When using the Dirac delta function approximation that is phase weighted"
                    "reciprocal times heaviside use positive weights! Abort..."));
      AssertThrow(std::abs(w_1g - w_1h) > std::numeric_limits<number>::epsilon(),
                  dealii::ExcMessage(
                    "When using the Dirac delta function approximation that is phase weighted"
                    "reciprocal times heaviside the phase weights must differ! Abort..."));
      AssertThrow(w_2g > 0.0,
                  dealii::ExcMessage(
                    "When using the Dirac delta function approximation that is phase weighted"
                    "reciprocal times heaviside use positive weights! Abort..."));
      AssertThrow(std::abs(w_2g - w_2h) > std::numeric_limits<number>::epsilon(),
                  dealii::ExcMessage(
                    "When using the Dirac delta function approximation that is phase weighted"
                    "reciprocal times heaviside the phase weights must differ! Abort..."));
    }
 
    inline number
    compute_weight(const number level_set_heaviside) const override
    {
      return compute_weight_internal(level_set_heaviside);
    }
 
    inline dealii::VectorizedArray<number>
    compute_weight(const dealii::VectorizedArray<number> &level_set_heaviside) const override
    {
      return compute_weight_internal(level_set_heaviside);
    }
 
  private:
    template <typename value_type>
    inline value_type
    compute_weight_internal(const value_type &level_set_heaviside) const
    {
      const auto w_1 = Tools::interpolate_reciprocal(level_set_heaviside, w_1g, w_1h);
      const auto w_2 = Tools::interpolate(level_set_heaviside, w_2g, w_2h);
      return w_1 * w_2 * correction_factor;
    }
 
    const number w_1g;
    const number w_1h;
    const number w_2g;
    const number w_2h;
 
    /*
     * correction factor
     *
     *           /      w_1g w_1h ln( w_1g / w_1h )
     * c_corr = | w_2g -----------------------------
     *          \              w_1g - w_1h
     *
     *                             1 / w_1g ( ln( w_1h / w_1g ) - 1 ) + 1 / w_1h  \-1
     *          + (w_2h - w_2g) * ----------------------------------------------- |
     *                                         ( 1/w_1h - 1/w_1g )²              /
     *
     */
    const number correction_factor;
  };
 
  template <typename number>
  std::unique_ptr<DeltaApproximationBase<number>>
  create_phase_weighted_delta_approximation(const DeltaApproximationPhaseWeightedData<number> &data)
  {
    switch (data.type)
      {
        case DiracDeltaFunctionApproximationType::norm_of_indicator_gradient:
          return nullptr;
        case DiracDeltaFunctionApproximationType::heaviside_phase_weighted:
          return std::make_unique<DeltaApproximationHeavisidePhaseWeighted<number>>(data);
        case DiracDeltaFunctionApproximationType::heavy_phase_only:
          return std::make_unique<DeltaApproximationHeavyPhaseOnly<number>>();
        case DiracDeltaFunctionApproximationType::heaviside_times_heaviside_phase_weighted:
          return std::make_unique<DeltaApproximationHeavisideTimesHeavisidePhaseWeighted<number>>(
            data);
        case DiracDeltaFunctionApproximationType::reciprocal_phase_weighted:
          return std::make_unique<DeltaApproximationReciprocalPhaseWeighted<number>>(data);
        case DiracDeltaFunctionApproximationType::reciprocal_times_heaviside_phase_weighted:
          return std::make_unique<DeltaApproximationReciprocalTimesHeavisidePhaseWeighted<number>>(
            data);
        default:
          Assert(false, dealii::ExcNotImplemented());
      }
    return nullptr;
  }
} // namespace MeltPoolDG::LevelSet