dg_generic_convection_diffusion_worker.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/tensor.h>
#include <deal.II/base/vectorization.h>
 
#include <meltpooldg/utilities/dealii_tensor.hpp>
 
namespace MeltPoolDG::Utils
{
 
  template <typename T>
  concept ConvectionKernelType = requires(const T &kernel) {
    typename T::ValueType;
    typename T::PhysicalFluxType;
 
    {
      kernel.flux(std::declval<typename T::ValueType>())
    } -> std::same_as<typename T::PhysicalFluxType>;
 
    {
      kernel.local_maximum_wave_speed(std::declval<typename T::ValueType>(),
                                      std::declval<typename T::ValueType>())
    };
  };
 
  template <int dim,
            typename number,
            ConvectionKernelType Kernel,
            typename VectorizedArrayType = dealii::VectorizedArray<number>>
  struct DGConvectionOperator
  {
    using ValueType = typename Kernel::ValueType;
    using FluxType  = typename Kernel::PhysicalFluxType;
 
    struct FaceFlux
    {
      ValueType inner_face_value;
      ValueType outer_face_value;
    };
 
    [[nodiscard]] static FluxType
    cell(const ValueType &u, const Kernel &kernel)
    {
      return kernel.flux(u);
    }
 
    [[nodiscard]] static FaceFlux
    face(const ValueType                                   &u_m,
         const ValueType                                   &u_p,
         const dealii::Tensor<1, dim, VectorizedArrayType> &normal,
         const Kernel                                      &kernel)
    {
      const FluxType flux_m = kernel.flux(u_m);
      const FluxType flux_p = kernel.flux(u_p);
 
      const auto local_maximum_wave_speed = kernel.local_maximum_wave_speed(u_m, u_p);
 
      const auto lax_friedrichs_flux = contract_average_tensor_with_vector(flux_m, flux_p, normal) +
                                       0.5 * local_maximum_wave_speed * (u_m - u_p);
 
      return FaceFlux{.inner_face_value = -lax_friedrichs_flux,
                      .outer_face_value = lax_friedrichs_flux};
    }
  };
 
  template <typename T>
  concept DiffusionKernelType = requires(const T &kernel) {
    typename T::ValueType;
    typename T::GradientType;
 
    {
      kernel.flux(std::declval<typename T::ValueType>(), std::declval<typename T::GradientType>())
    } -> std::same_as<typename T::GradientType>;
  };
 
  template <int dim,
            typename number,
            DiffusionKernelType Kernel,
            typename VectorizedArrayType = dealii::VectorizedArray<number>>
  struct DGDiffusionOperator
  {
    using ValueType    = typename Kernel::ValueType;
    using GradientType = typename Kernel::GradientType;
 
    struct FaceFlux
    {
      ValueType    inner_face_value;
      ValueType    outer_face_value;
      GradientType inner_face_gradient;
      GradientType outer_face_gradient;
    };
 
    [[nodiscard]] static GradientType
    cell(const ValueType &u, const GradientType &grad_u, const Kernel &kernel)
    {
      return -1. * kernel.flux(u, grad_u);
    }
 
    [[nodiscard]] static FaceFlux
    face(const ValueType                                   &u_m,
         const ValueType                                   &u_p,
         const GradientType                                &grad_u_m,
         const GradientType                                &grad_u_p,
         const dealii::Tensor<1, dim, VectorizedArrayType> &normal,
         const VectorizedArrayType                          penalty_parameter,
         const Kernel                                      &kernel)
    {
      GradientType flux_m = kernel.flux(u_m, grad_u_m);
      GradientType flux_p = kernel.flux(u_p, grad_u_p);
 
      const auto value = contract_average_tensor_with_vector(flux_m, flux_p, normal) -
                         penalty_parameter * (u_m - u_p);
 
 
      GradientType jump_u = jump(u_m, u_p, normal);
 
      // use jumps instead of gradients for evaluating the diffusive flux
      const GradientType penalty_flux_m = 0.5 * kernel.flux(u_m, jump_u);
      const GradientType penalty_flux_p = 0.5 * kernel.flux(u_p, jump_u);
 
      return FaceFlux{.inner_face_value    = value,
                      .outer_face_value    = -value,
                      .inner_face_gradient = penalty_flux_m,
                      .outer_face_gradient = penalty_flux_p};
    }
  };
 
  template <int dim,
            typename number,
            ConvectionKernelType ConvectiveKernel,
            DiffusionKernelType  DiffusiveKernel,
            typename VectorizedArrayType = dealii::VectorizedArray<number>>
  struct DGConvectionDiffusionOperator
  {
    using ConvectionOperator =
      DGConvectionOperator<dim, number, ConvectiveKernel, VectorizedArrayType>;
    using DiffusionOperator =
      DGDiffusionOperator<dim, number, DiffusiveKernel, VectorizedArrayType>;
 
    using ValueType    = typename DiffusionOperator::ValueType;
    using GradientType = typename DiffusionOperator::GradientType;
 
    struct FaceFlux
    {
      ValueType    inner_face_value;
      ValueType    outer_face_value;
      GradientType inner_face_gradient;
      GradientType outer_face_gradient;
    };
 
    [[nodiscard]] static GradientType
    cell(const ValueType        &u,
         const GradientType     &grad_u,
         const ConvectiveKernel &convective_kernel,
         const DiffusiveKernel  &diffusion_kernel)
    {
      GradientType flux = ConvectionOperator::cell(u, convective_kernel);
      flux += DiffusionOperator::cell(u, grad_u, diffusion_kernel);
      return flux;
    }
 
    [[nodiscard]] static FaceFlux
    face(const ValueType                                   &u_m,
         const ValueType                                   &u_p,
         const GradientType                                &grad_u_m,
         const GradientType                                &grad_u_p,
         const dealii::Tensor<1, dim, VectorizedArrayType> &normal,
         const VectorizedArrayType                          penalty_parameter,
         const ConvectiveKernel                            &convective_kernel,
         const DiffusiveKernel                             &diffusion_kernel)
    {
      const auto convection_face_value =
        ConvectionOperator::face(u_m, u_p, normal, convective_kernel);
      const auto diffusion_face_value = DiffusionOperator::face(
        u_m, u_p, grad_u_m, grad_u_p, normal, penalty_parameter, diffusion_kernel);
 
      return FaceFlux{
        .inner_face_value =
          convection_face_value.inner_face_value + diffusion_face_value.inner_face_value,
        .outer_face_value =
          convection_face_value.outer_face_value + diffusion_face_value.outer_face_value,
        .inner_face_gradient = diffusion_face_value.inner_face_gradient,
        .outer_face_gradient = diffusion_face_value.outer_face_gradient,
      };
    }
  };
 
} // namespace MeltPoolDG::Utils