multiphase_operator.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/fe/fe_system.h>
 
#include <meltpooldg/compressible_flow/convective_kernels.hpp>
#include <meltpooldg/compressible_flow/data_types.hpp>
#include <meltpooldg/compressible_flow/kernels.hpp>
#include <meltpooldg/compressible_flow/multiphase_level_set_advection.hpp>
#include <meltpooldg/cut/util.hpp>
#include <meltpooldg/flow/darcy_damping_model.hpp>
#include <meltpooldg/phase_change/evaporation_model_knight.hpp>
#include <meltpooldg/utilities/better_enum.hpp>
#include <meltpooldg/utilities/dg_generic_convection_diffusion_worker.hpp>
#include <meltpooldg/utilities/vector_tools.hpp>
 
namespace MeltPoolDG::Multiphase
{
  BETTER_ENUM(Phase, int, liquid, gas);
 
  template <int dim, typename number, bool is_viscous_gas = true, bool is_viscous_liquid = true>
  class CompressibleMultiphaseOperator
  {
  public:
    using VectorType            = dealii::LinearAlgebra::distributed::Vector<number>;
    using MappingInfoType       = CutUtil::MappingInfoType<dim, number>;
    using MappingInfoVectorType = CutUtil::MappingInfoVectorType<dim, number>;
 
    using ConservedVariablesType = CompressibleFlow::ConservedVariablesType<dim, number>;
    using ConservedVariablesGradType =
      CompressibleFlow::ConservedVariablesGradientType<dim, number>;
 
    using FlowFluxType   = CompressibleFlow::FluxType<dim, number>;
    using FlowSourceType = CompressibleFlow::SourceType<dim, number>;
 
    using ConvectiveKernel = CompressibleFlow::
      ConvectiveFlux<dim, number, ConservedVariablesType, ConservedVariablesGradType>;
    using DiffusiveKernel = CompressibleFlow::
      DiffusiveFlux<dim, number, ConservedVariablesType, ConservedVariablesGradType, FlowFluxType>;
 
    using ConvectionDiffusionOperator =
      Utils::DGConvectionDiffusionOperator<dim, number, ConvectiveKernel, DiffusiveKernel>;
 
    using ConvectionOperator = Utils::DGConvectionOperator<dim, number, ConvectiveKernel>;
 
    using DofValueView = CompressibleFlow::DofValueView<dim, ConservedVariablesType>;
    using DofPrimitiveValueView =
      CompressibleFlow::DofPrimitiveValueView<dim, ConservedVariablesType>;
    using DofPrimitiveStateView =
      CompressibleFlow::DofPrimitiveStateView<dim, number, const ConservedVariablesType>;
    using DofStateView = CompressibleFlow::DofStateView<dim, number, const ConservedVariablesType>;
    using DofValueAndGradientStateView =
      CompressibleFlow::DofValueAndGradientStateView<dim,
                                                     number,
                                                     const ConservedVariablesType,
                                                     const ConservedVariablesGradType>;
 
    template <int n_components = CompressibleFlow::n_conserved_variables<dim>>
    using DomainEval = FECellIntegrator<dim, n_components, number>;
 
    template <int n_components = CompressibleFlow::n_conserved_variables<dim>>
    using FaceEval = FEFaceIntegrator<dim, n_components, number>;
 
    template <int n_components = CompressibleFlow::n_conserved_variables<dim>>
    using DomainPointEval =
      dealii::FEPointEvaluation<n_components, dim, dim, dealii::VectorizedArray<number>>;
 
    template <int n_components = CompressibleFlow::n_conserved_variables<dim>>
    using FacePointEval =
      dealii::FEFacePointEvaluation<n_components, dim, dim, dealii::VectorizedArray<number>>;
 
    explicit CompressibleMultiphaseOperator(
      CompressibleFlow::MultiphaseOperationScratchData<dim, number> &multiphase_scratch_data,
      const MappingInfoType                                         &mapping_info_interface_in,
      const MappingInfoVectorType                                   &mapping_info_cells_in,
      const MappingInfoVectorType                                   &mapping_info_faces_in);
 
    void
    local_apply_cell_rhs(const dealii::MatrixFree<dim, number>       &matrix_free,
                         VectorType                                  &dst,
                         const VectorType                            &src,
                         const std::pair<unsigned int, unsigned int> &cell_range) const;
 
    void
    local_apply_face_rhs(const dealii::MatrixFree<dim, number>       &matrix_free,
                         VectorType                                  &dst,
                         const VectorType                            &src,
                         const std::pair<unsigned int, unsigned int> &face_range) const;
 
    void
    local_apply_boundary_face_rhs(const dealii::MatrixFree<dim, number>       &matrix_free,
                                  VectorType                                  &dst,
                                  const VectorType                            &src,
                                  const std::pair<unsigned int, unsigned int> &face_range) const;
 
    void
    local_apply_cell_lhs(const dealii::MatrixFree<dim, number>       &matrix_free,
                         VectorType                                  &dst,
                         const VectorType                            &src,
                         const std::pair<unsigned int, unsigned int> &cell_range) const;
 
    void
    local_apply_face_lhs(const dealii::MatrixFree<dim, number>       &matrix_free,
                         VectorType                                  &dst,
                         const VectorType                            &src,
                         const std::pair<unsigned int, unsigned int> &face_range) const;
 
    void
    local_apply_boundary_face_lhs(const dealii::MatrixFree<dim, number>       &matrix_free,
                                  VectorType                                  &dst,
                                  const VectorType                            &src,
                                  const std::pair<unsigned int, unsigned int> &face_range) const;
 
    void
    create_rhs(const number     &time,
               const number     &time_step_in,
               VectorType       &dst,
               const VectorType &src) const;
 
    void
    vmult(VectorType &dst, const VectorType &src) const;
 
    void
    set_current_time(const number &current_time_in)
    {
      current_time = current_time_in;
    };
 
    void
    add_external_force(
      std::shared_ptr<CompressibleFlow::ExternalFlowForce<dim, number>> external_force);
 
    template <typename EvaluatorType>
    inline void
    ghost_penalty_face_integral(EvaluatorType     &eval_m,
                                EvaluatorType     &eval_p,
                                const unsigned int q,
                                const number       cell_side_length,
                                const number       cell_side_length_pow_3,
                                const number       cell_side_length_pow_5) const;
 
    void
    compute_inverse_diagonal_from_matrixfree(VectorType &diagonal) const;
 
    void
    compute_system_matrix_from_matrixfree(
      dealii::TrilinosWrappers::SparseMatrix &system_matrix) const;
 
  private:
    CompressibleFlow::MultiphaseOperationScratchData<dim, number> &multiphase_scratch_data;
 
    const MappingInfoType &mapping_info_interface;
 
    const MappingInfoVectorType &mapping_info_cells;
 
    const MappingInfoVectorType &mapping_info_faces;
 
    dealii::FESystem<dim> fe_point_temp;
 
    const unsigned int n_dofs_per_cell;
 
    number visc_ave_weight_phase_liquid;
    number visc_ave_weight_phase_gas;
 
    std::vector<std::shared_ptr<CompressibleFlow::ExternalFlowForce<dim, number>>> external_forces;
 
    mutable number time_step = 0.;
 
    mutable number inv_time_step = 0.;
 
    mutable number current_time = 0.;
 
    mutable LevelSetAdvection<dim, number> level_set_advection_operator;
 
    std::unique_ptr<Evaporation::EvaporationModelKnight<number>> evaporation_model_knight;
 
    Flow::DarcyDampingModel<number> darcy_damping_model;
 
    static constexpr std::array<unsigned int, 2> first_component = {{
      0 /*liquid phase*/,
      CompressibleFlow::n_conserved_variables<dim> /*gas phase*/
    }};
 
    FECellIntegrator<dim, dim + 2, number>
    create_cell_integrator(const CutUtil::CellCategory category,
                           const unsigned int          offset = 0) const
    {
      return FECellIntegrator<dim, dim + 2, number>(
        multiphase_scratch_data.scratch_data.get_matrix_free(),
        multiphase_scratch_data.dof_idx,
        multiphase_scratch_data.quad_idx,
        offset,
        category);
    }
 
    FEFaceIntegrator<dim, dim + 2, number>
    create_face_integrator(const bool                  is_inner_face,
                           const CutUtil::CellCategory category,
                           const unsigned int          offset = 0) const
    {
      return FEFaceIntegrator<dim, dim + 2, number>(
        multiphase_scratch_data.scratch_data.get_matrix_free(),
        is_inner_face,
        multiphase_scratch_data.dof_idx,
        multiphase_scratch_data.quad_idx,
        offset,
        category);
    };
 
    std::pair<FEFaceIntegrator<dim, dim + 2, number>, FEFaceIntegrator<dim, dim + 2, number>>
    create_face_integrators(const CutUtil::CellCategory category,
                            const unsigned int          offset = 0) const
    {
      return {create_face_integrator(true, category, offset),
              create_face_integrator(false, category, offset)};
    };
 
    void
    internal_compute_diagonal_or_system_matrix(
      VectorType                             &diagonal,
      dealii::TrilinosWrappers::SparseMatrix &system_matrix,
      const bool                              do_diagonal) const;
  };
} // namespace MeltPoolDG::Multiphase