utils.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <deal.II/base/aligned_vector.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/vectorization.h>
 
#include <deal.II/grid/reference_cell.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
 
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/la_parallel_vector.h>
 
#include "meltpooldg/compressible_flow/material.hpp"
#include <meltpooldg/compressible_flow/data_types.hpp>
#include <meltpooldg/core/scratch_data.hpp>
#include <meltpooldg/utilities/better_enum.hpp>
#include <meltpooldg/utilities/matrix_free_util.hpp>
 
namespace MeltPoolDG::CompressibleFlow
{
  // Forward declaration
  template <int dim, typename number>
  class Material;
 
  BETTER_ENUM(Idx1D, char, density, momentum_x, energy);
  BETTER_ENUM(Idx2D, char, density, momentum_x, momentum_y, energy);
  BETTER_ENUM(Idx3D, char, density, momentum_x, momentum_y, momentum_z, energy);
 
  template <int dim, typename number>
  inline DEAL_II_ALWAYS_INLINE //
    dealii::Tensor<1, dim, dealii::VectorizedArray<number>>
    calculate_velocity(const ConservedVariablesType<dim, number> &conserved_variables);
 
  template <int dim, typename number>
  inline DEAL_II_ALWAYS_INLINE //
    dealii::Tensor<2, dim, dealii::VectorizedArray<number>>
    calculate_grad_velocity(
      const ConservedVariablesType<dim, number>         &conserved_variables,
      const ConservedVariablesGradientType<dim, number> &grad_conserved_variables);
 
  template <int dim, typename Number>
  void
  calculate_penalty_parameter(
    dealii::AlignedVector<dealii::VectorizedArray<Number>> &array_penalty_parameter,
    const dealii::MatrixFree<dim, Number>                  &matrix_free,
    const std::string                                      &domain_representation_type,
    const unsigned int                                      dof_index      = 0,
    const Number                                            scaling_factor = 1.0);
 
  template <int dim, typename number>
  void
  update_primitive_variables_solution(
    dealii::LinearAlgebra::distributed::Vector<number>       &solution_primitive_variables,
    const dealii::LinearAlgebra::distributed::Vector<number> &solution,
    const ScratchData<dim, dim, number>                      &scratch_data,
    const unsigned int                                        dof_idx,
    const unsigned int                                        quad_idx,
    const Material<dim, number>                              *material_liquid,
    const Material<dim, number>                              *material_gas = nullptr);
 
 
  template <typename number, int n_species>
  inline bool
  is_viscous_flow(const MaterialPhaseData<number> &material_data);
 
  template <int dim, typename number, int n_species = 1>
  struct ExternalFlowForce
  {
    virtual ~ExternalFlowForce() = default;
 
    virtual ConservedVariablesType<dim, number, n_species>
    value(number                                                     time_step_size,
          unsigned int                                               cell_batch_id,
          const dealii::Point<dim, dealii::VectorizedArray<number>> &points,
          const ConservedVariablesType<dim, number, n_species>      &w) = 0;
  };
 
 
  template <int dim, typename number, int n_species = 1>
  struct ExternalFlowForceJacobian
  {
    virtual ~ExternalFlowForceJacobian() = default;
 
    virtual ConservedVariablesType<dim, number, n_species>
    value(number                                                     time_step_size,
          unsigned int                                               cell_batch_id,
          const dealii::Point<dim, dealii::VectorizedArray<number>> &points,
          const ConservedVariablesType<dim, number, n_species>      &w,
          const ConservedVariablesType<dim, number, n_species>      &delta_w) = 0;
  };
 
  /********************************************************************************************
   * Inlined function definitions
   * *************************************************************************************+****/
  template <int dim, typename number>
  inline DEAL_II_ALWAYS_INLINE //
    dealii::Tensor<1, dim, dealii::VectorizedArray<number>>
    calculate_velocity(const ConservedVariablesType<dim, number> &conserved_variables)
  {
    const dealii::VectorizedArray<number> inverse_density =
      dealii::VectorizedArray<number>(1.) / conserved_variables[0];
 
    dealii::Tensor<1, dim, dealii::VectorizedArray<number>> velocity;
    for (unsigned int d = 0; d < dim; ++d)
      velocity[d] = conserved_variables[1 + d] * inverse_density;
 
    return velocity;
  }
 
  template <int dim, typename number>
  inline DEAL_II_ALWAYS_INLINE //
    dealii::Tensor<2, dim, dealii::VectorizedArray<number>>
    calculate_grad_velocity(
      const ConservedVariablesType<dim, number>         &conserved_variables,
      const ConservedVariablesGradientType<dim, number> &grad_conserved_variables)
  {
    const dealii::VectorizedArray<number> inverse_density =
      dealii::VectorizedArray<number>(1.) / conserved_variables[0];
    const dealii::Tensor<1, dim, dealii::VectorizedArray<number>> velocity =
      calculate_velocity<dim, number>(conserved_variables);
 
    dealii::Tensor<1, dim, dealii::VectorizedArray<number>> grad_rho;
    for (unsigned int d = 0; d < dim; ++d)
      grad_rho[d] = grad_conserved_variables[0][d];
 
    dealii::Tensor<2, dim, dealii::VectorizedArray<number>> grad_rho_velocity;
    for (unsigned int d = 0; d < dim; ++d)
      for (unsigned int e = 0; e < dim; ++e)
        grad_rho_velocity[d][e] = grad_conserved_variables[1 + d][e];
 
    dealii::Tensor<2, dim, dealii::VectorizedArray<number>> grad_velocity;
    for (unsigned int d = 0; d < dim; ++d)
      for (unsigned int e = 0; e < dim; ++e)
        grad_velocity[d][e] =
          inverse_density * (grad_rho_velocity[d][e] - velocity[d] * grad_rho[e]);
 
    return grad_velocity;
  }
 
 
  template <typename number, int n_species>
  inline bool
  is_viscous_flow(const MaterialPhaseData<number> &material_data)
  {
    for (unsigned int species = 0; species < n_species; ++species)
      if (material_data.species_data[species].dynamic_viscosity > 0.0)
        return true;
    return false;
  }
 
  template <int dim, typename Number>
  void
  calculate_penalty_parameter(
    dealii::AlignedVector<dealii::VectorizedArray<Number>> &array_penalty_parameter,
    const dealii::MatrixFree<dim, Number>                  &matrix_free,
    const std::string                                      &domain_representation_type,
    const unsigned int                                      dof_index,
    const Number                                            scaling_factor)
  {
    const unsigned int n_cells = matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();
    array_penalty_parameter.resize(n_cells);
 
    dealii::Mapping<dim> const       &mapping = *matrix_free.get_mapping_info().mapping;
    dealii::FiniteElement<dim> const &fe      = matrix_free.get_dof_handler(dof_index).get_fe();
    unsigned int const                degree  = fe.degree;
 
    // use penalty factor for hypercube elements according to K. Hillewaert, Development of the
    // discontinuous Galerkin method for high-resolution, large scale CFD and acoustics in
    // industrial geometries, PhD thesis, Univ. de Louvain, 2013.
    const Number fac = scaling_factor * (degree + 1.0) * (degree + 1.0);
 
    const std::vector<dealii::ReferenceCell<dim>> reference_cells =
      matrix_free.get_dof_handler(dof_index).get_triangulation().get_reference_cells();
    AssertThrow(reference_cells.size() == 1, dealii::ExcMessage("No mixed meshes allowed."));
 
 
    const dealii::Quadrature<dim> quadrature =
      reference_cells[0].get_gauss_type_quadrature(degree + 1);
    dealii::FEValues<dim> fe_values(mapping, fe, quadrature, dealii::update_JxW_values);
 
    const dealii::Quadrature<dim - 1> face_quadrature =
      reference_cells[0].face_reference_cell(0).get_gauss_type_quadrature(degree + 1);
    dealii::FEFaceValues<dim> fe_face_values(mapping,
                                             fe,
                                             face_quadrature,
                                             dealii::update_JxW_values);
 
    if (domain_representation_type == "fitted")
      {
        for (unsigned int i = 0; i < n_cells; ++i)
          {
            for (unsigned int v = 0; v < matrix_free.n_active_entries_per_cell_batch(i); ++v)
              {
                typename dealii::DoFHandler<dim>::cell_iterator cell =
                  matrix_free.get_cell_iterator(i, v, dof_index);
                fe_values.reinit(cell);
 
                // calculate cell volume
                Number volume = 0;
                for (unsigned int q = 0; q < quadrature.size(); ++q)
                  {
                    volume += fe_values.JxW(q);
                  }
 
                // calculate surface area
                Number surface_area = 0;
                for (unsigned int const f : cell->face_indices())
                  {
                    fe_face_values.reinit(cell, f);
                    Number const factor =
                      (cell->at_boundary(f) and not(cell->has_periodic_neighbor(f))) ? 1. : 0.5;
                    for (unsigned int q = 0; q < face_quadrature.size(); ++q)
                      {
                        surface_area += fe_face_values.JxW(q) * factor;
                      }
                  }
 
                array_penalty_parameter[i][v] = surface_area / volume * fac;
              }
          }
      }
    else if (domain_representation_type == "cut")
      {
        for (unsigned int i = 0; i < n_cells; ++i)
          {
            for (unsigned int v = 0; v < matrix_free.n_active_entries_per_cell_batch(i); ++v)
              {
                typename dealii::DoFHandler<dim>::cell_iterator cell =
                  matrix_free.get_cell_iterator(i, v, dof_index);
 
                // simplified computation for hypercube elements
                array_penalty_parameter[i][v] = fac / cell->minimum_vertex_distance();
              }
          }
      }
    else
      AssertThrow(false,
                  dealii::ExcMessage("The domain representation type '" +
                                     domain_representation_type + "' is not supported."));
  }
 
  template <int dim, typename number>
  void
  update_primitive_variables_solution(
    dealii::LinearAlgebra::distributed::Vector<number>       &solution_primitive_variables,
    const dealii::LinearAlgebra::distributed::Vector<number> &solution,
    const ScratchData<dim, dim, number>                      &scratch_data,
    const unsigned int                                        dof_idx,
    const unsigned int                                        quad_idx,
    const Material<dim, number>                              *material_liquid,
    const Material<dim, number>                              *material_gas)
  {
    const dealii::MatrixFree<dim, number, dealii::VectorizedArray<number>> &matrix_free =
      scratch_data.get_matrix_free();
    unsigned int n_support_points_per_cell = scratch_data.get_n_dofs_per_cell(dof_idx) / (dim + 2);
 
    auto process_cell = [&](CutUtil::CellCategory        category,
                            unsigned int                 first_component,
                            const Material<dim, number> *material,
                            const unsigned int           cell_batch) {
      if (!material)
        return;
      auto eval = FECellIntegrator<dim, dim + 2, number>(
        matrix_free, dof_idx, quad_idx, first_component, category);
      eval.reinit(cell_batch);
      eval.read_dof_values(solution);
 
      for (unsigned int i = 0; i < n_support_points_per_cell; ++i)
        {
          const auto &u_cons = eval.get_dof_value(i);
          auto u_prim = material->eos_utils->convert_conservative_into_primitive_variables(u_cons);
          eval.submit_dof_value(u_prim, i);
        }
 
      eval.set_dof_values(solution_primitive_variables);
    };
 
    if (!matrix_free.get_dof_info(dof_idx).cell_active_fe_index.empty())
      {
        // using hp::FECollection (cutFEM/DG)
        for (unsigned int cell_batch = 0; cell_batch < matrix_free.n_cell_batches(); ++cell_batch)
          {
            const auto cell_category = matrix_free.get_cell_category(cell_batch);
 
            if (cell_category == CutUtil::CellCategory::liquid)
              {
                process_cell(CutUtil::CellCategory::liquid, 0, material_liquid, cell_batch);
              }
            else if (cell_category == CutUtil::CellCategory::gas)
              {
                process_cell(CutUtil::CellCategory::gas, dim + 2, material_gas, cell_batch);
              }
            else if (cell_category == CutUtil::CellCategory::intersected)
              {
                process_cell(CutUtil::CellCategory::intersected, 0, material_liquid, cell_batch);
                process_cell(CutUtil::CellCategory::intersected, dim + 2, material_gas, cell_batch);
              }
          }
      }
    else
      {
        // not using hp::FECollection (fitted mesh)
        for (unsigned int cell_batch = 0; cell_batch < matrix_free.n_cell_batches(); ++cell_batch)
          {
            process_cell(CutUtil::CellCategory::liquid, 0, material_liquid, cell_batch);
          }
      }
  }
 
  template <typename DofViewType, typename VectorizedArrayType>
  inline DEAL_II_ALWAYS_INLINE //
    VectorizedArrayType
    maximum_local_wave_speed(const DofViewType &u_m, const DofViewType &u_p)
  {
    const auto velocity_m = u_m.velocity();
    const auto velocity_p = u_p.velocity();
 
    const auto sound_speed_p = u_p.speed_of_sound();
    const auto sound_speed_m = u_m.speed_of_sound();
 
    const auto sound_speed_p2 = sound_speed_p * sound_speed_p;
    const auto sound_speed_m2 = sound_speed_m * sound_speed_m;
 
    const auto lambda = 0.5 * std::sqrt(std::max(velocity_p.norm_square() + sound_speed_p2,
                                                 velocity_m.norm_square() + sound_speed_m2));
    return lambda;
  }
} // namespace MeltPoolDG::CompressibleFlow