multiphase_level_set_advection.hpp

Go to the documentation of this file.

 
// TODO: remove/restructure, when level-set advection and reinitialization (and velocity field
// extrapolation) are incorporated
 
#pragma once
 
#include <deal.II/base/function.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/point.h>
#include <deal.II/base/utilities.h>
 
#include <deal.II/numerics/vector_tools_interpolate.h>
 
#include <numeric>
 
namespace MeltPoolDG::Multiphase
{
  template <int dim, typename number>
  class LevelSetOscillatingWaterColumn : public dealii::Function<dim, number>
  {
  public:
    explicit LevelSetOscillatingWaterColumn()
      : dealii::Function<dim, number>(1)
    {}
 
    void
    update_interface_positions(const number              &time_step,
                               const std::vector<number> &interface_velocity)
    {
      for (unsigned int i = 0; i < 2; ++i)
        interface_coordinates[i] += time_step * interface_velocity[i];
    }
 
    number
    value(const dealii::Point<dim, number> &p, const unsigned int /*component*/) const final
    {
      return std::min(p[0] - interface_coordinates[0], -(p[0] - interface_coordinates[1]));
    }
 
  private:
    // Note that the initial interface coordinates are aligned to the data in the case file
    // "oscillating_water_column.hpp".
    std::vector<number> interface_coordinates = {0.100156487, 0.800156487};
  };
 
  template <int dim, typename number>
  class LevelSetAdvection
  {
  public:
    LevelSetAdvection();
 
    static LevelSetOscillatingWaterColumn<dim, number> level_set_oscillating_water_column;
 
    void
    move_level_set(auto              &level_set,
                   const number      &time_step,
                   const std::string &case_name,
                   const auto        &scratch_data,
                   const int          level_set_dof_idx)
    {
      if (case_name == "oscillating_water_column")
        {
          // two phase boundaries
          level_set_oscillating_water_column.update_interface_positions(time_step,
                                                                        interface_velocity);
          level_set.zero_out_ghost_values();
          dealii::VectorTools::interpolate(scratch_data->get_mapping(),
                                           scratch_data->get_dof_handler(level_set_dof_idx),
                                           level_set_oscillating_water_column,
                                           level_set);
          level_set.update_ghost_values();
        }
      else if (case_name == "two_phase" or case_name == "shell")
        {
          // one phase boundary
 
          // TODO: This is a current work-around, since the level-set advection is not yet fully
          // implemented. The mean value of the interface velocity is used to advect the level-set
          // field. In the future, a full level-set advection should be implemented, which uses the
          // full velocity field to advect the level-set field.
          const number local_sum =
            std::accumulate(interface_velocity.begin(), interface_velocity.end(), 0.);
 
          const unsigned int local_size = interface_velocity.size();
 
          const number global_sum =
            dealii::Utilities::MPI::sum(local_sum, scratch_data->get_mpi_comm(level_set_dof_idx));
 
          const unsigned int global_size =
            dealii::Utilities::MPI::sum(local_size, scratch_data->get_mpi_comm(level_set_dof_idx));
 
          const number velocity = global_sum / global_size;
 
          for (unsigned int i = 0; i < level_set.locally_owned_size(); ++i)
            level_set.local_element(i) += time_step * velocity;
 
          level_set.update_ghost_values();
        }
      else
        AssertThrow(false,
                    dealii::ExcMessage(
                      "Analytical function for level-set advection is "
                      "only supported for the cases 'oscillating_water_column', 'two_phase', "
                      "and 'shell'."));
    }
 
    void
    set_interface_velocity(const number &interface_velocity_in)
    {
      interface_velocity.push_back(interface_velocity_in);
    }
 
    void
    clear_interface_velocity()
    {
      interface_velocity.clear();
    }
 
  private:
    static std::vector<number> interface_velocity;
  };
 
  template <int dim, typename number>
  LevelSetAdvection<dim, number>::LevelSetAdvection() = default;
 
  template <int dim, typename number>
  std::vector<number> LevelSetAdvection<dim, number>::interface_velocity;
 
  template <int dim, typename number>
  LevelSetOscillatingWaterColumn<dim, number>
    LevelSetAdvection<dim, number>::level_set_oscillating_water_column;
} // namespace MeltPoolDG::Multiphase