incompressible_newtonian_evaporation_material.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/tensor.h>
#include <deal.II/base/vectorization.h>
 
#include <deal.II/lac/la_parallel_block_vector.h>
#include <deal.II/lac/la_parallel_vector.h>
 
#include <deal.II/matrix_free/evaluation_flags.h>
 
#include <meltpooldg/core/scratch_data.hpp>
#include <meltpooldg/flow/incompressible_flow_material_base.hpp>
#include <meltpooldg/utilities/dealii_tensor.hpp>
#include <meltpooldg/utilities/fe_integrator.hpp>
 
#include <functional>
 
 
namespace MeltPoolDG::Evaporation
{
  template <int dim, typename number>
  class IncompressibleNewtonianFluidEvaporationMaterial
    : public Flow::IncompressibleMaterialBase<dim, number>
  {
  private:
    using BlockVectorType = dealii::LinearAlgebra::distributed::BlockVector<number>;
 
  public:
    IncompressibleNewtonianFluidEvaporationMaterial(
      const ScratchData<dim, dim, number>                                    &scratch_data,
      const std::function<const dealii::VectorizedArray<number>
                            &(const unsigned int cell, const unsigned int q)> get_viscosity,
      const BlockVectorType                                                  &normal_vector,
      const dealii::LinearAlgebra::distributed::Vector<number>               &heaviside,
      const unsigned int                                                      normal_dof_idx,
      const unsigned int ls_hanging_nodes_dof_idx,
      const unsigned int velocity_quad_idx)
      : scratch_data(scratch_data)
      , get_viscosity(get_viscosity)
      , normal_vector(normal_vector)
      , heaviside(heaviside)
      , normal_dof_idx(normal_dof_idx)
      , ls_hanging_nodes_dof_idx(ls_hanging_nodes_dof_idx)
      , velocity_quad_idx(velocity_quad_idx)
      , normal_vals(scratch_data.get_matrix_free(), normal_dof_idx, velocity_quad_idx)
      , ls_vals(scratch_data.get_matrix_free(), ls_hanging_nodes_dof_idx, velocity_quad_idx)
    {}
 
    void
    reinit(const dealii::Tensor<2, dim, dealii::VectorizedArray<number>> &velocity_gradient,
           const unsigned int                                             cell_idx,
           const unsigned int                                             quad_idx) final
    {
      grad_u = velocity_gradient;
      div_u  = dealii::trace(velocity_gradient);
 
      viscosity = get_viscosity(cell_idx, quad_idx);
 
      // read normal vector values only for first quadrature point
      if (quad_idx == 0 || cell_idx != cell)
        {
          normal_vals.reinit(cell_idx);
          normal_vals.read_dof_values_plain(normal_vector);
          normal_vals.evaluate(dealii::EvaluationFlags::values);
 
          ls_vals.reinit(cell_idx);
          ls_vals.read_dof_values_plain(heaviside);
          ls_vals.evaluate(dealii::EvaluationFlags::values);
        }
 
      normal = normalize<dim>(normal_vals.get_value(quad_idx), 1e-10);
      hs     = ls_vals.get_value(quad_idx);
 
      cell = cell_idx;
    }
 
    dealii::Tensor<2, dim, dealii::VectorizedArray<number>>
    get_tau() final
    {
      const auto mask = dealii::compare_and_apply_mask<dealii::SIMDComparison::less_than>(
        hs,
        1.0,
        dealii::compare_and_apply_mask<dealii::SIMDComparison::greater_than>(hs, 0, 1.0, 0.0),
        0);
 
      return 2. * viscosity *
             (0.5 * (grad_u + transpose(grad_u)) - mask * div_u * outer_product(normal, normal));
    }
 
    dealii::Tensor<2, dim, dealii::VectorizedArray<number>>
    get_vmult_d_tau_d_grad_vel() final
    {
      // Since the velocity DoFs occur linear in the expression for tau, the vmult of the derivative
      // is the same operation, except with a different meaning of grad_u, i.e. shape functions
      // gradient times correction values for the velocity DoFs from the linear solver.
      return get_tau();
 
      // @note alternative, in case we want to convert this function to a real material tangent and
      // leave the multiplication with the velocity gradient to the external code.
      //
      // const auto identity  = Tensor<4,dim,VectorizedArray<number>>(identity_tensor<dim,
      // VectorizedArray<number>>()); const auto identity2 =
      // Tensor<2,dim,VectorizedArray<number>>(unit_symmetric_tensor<dim,
      // VectorizedArray<number>>());
 
      // const auto temp = 2. * viscosity * (identity - outer_product(identity2,
      // outer_product(normal, normal))); return double_contract<0,0,1,1>(temp, grad_u);
    }
 
    void
    update_ghost_values() final
    {
      normal_update_ghosts = !normal_vector.has_ghost_elements();
      if (normal_update_ghosts)
        normal_vector.update_ghost_values();
 
      ls_update_ghosts = !heaviside.has_ghost_elements();
      if (ls_update_ghosts)
        heaviside.update_ghost_values();
    }
 
    void
    zero_out_ghost_values() final
    {
      if (normal_update_ghosts)
        normal_vector.zero_out_ghost_values();
 
      if (ls_update_ghosts)
        heaviside.zero_out_ghost_values();
    }
 
  private:
    const ScratchData<dim, dim, number> &scratch_data;
    const std::function<const dealii::VectorizedArray<number> &(const unsigned int cell,
                                                                const unsigned int q)>
                                                              get_viscosity;
    const BlockVectorType                                    &normal_vector;
    const dealii::LinearAlgebra::distributed::Vector<number> &heaviside;
    const unsigned int                                        normal_dof_idx;
    const unsigned int                                        ls_hanging_nodes_dof_idx;
    const unsigned int                                        velocity_quad_idx;
    FECellIntegrator<dim, dim, number>                        normal_vals;
    FECellIntegrator<dim, 1, number>                          ls_vals;
 
    // temporary quadrature point values
    dealii::Tensor<2, dim, dealii::VectorizedArray<number>> grad_u;
    dealii::VectorizedArray<number>                         div_u;
    dealii::VectorizedArray<number>                         viscosity;
    dealii::Tensor<1, dim, dealii::VectorizedArray<number>> normal;
    dealii::VectorizedArray<number>                         hs;
    unsigned int                                            cell;
 
    mutable bool ls_update_ghosts     = true;
    mutable bool normal_update_ghosts = true;
  };
} // namespace MeltPoolDG::Evaporation