olb::FreeSurfaceMassExcessPostProcessor3D< T, DESCRIPTOR > Class Template Reference

Free Surface Processor 6 Calculates mass excess from the cell type conversions and distributes them to neighbouring interface cells Keeps mass local if no neighbour exists until an interface reappears at this position. More...

#include <freeSurfacePostProcessor3D.h>

Collaboration diagram for olb::FreeSurfaceMassExcessPostProcessor3D< T, DESCRIPTOR >:

Public Member Functions
int	getPriority () const
template<typename CELL >
void	apply (CELL &cell) any_platform
	Free Surface Processor 6: Mass Excess.

Static Public Attributes
static constexpr OperatorScope	scope = OperatorScope::PerCell

Detailed Description

template<typename T, typename DESCRIPTOR>
class olb::FreeSurfaceMassExcessPostProcessor3D< T, DESCRIPTOR >

Free Surface Processor 6 Calculates mass excess from the cell type conversions and distributes them to neighbouring interface cells Keeps mass local if no neighbour exists until an interface reappears at this position.

Definition at line 128 of file freeSurfacePostProcessor3D.h.

Member Function Documentation

apply()

template<typename T , typename DESCRIPTOR >

template<typename CELL >

void olb::FreeSurfaceMassExcessPostProcessor3D< T, DESCRIPTOR >::apply

(

CELL &

cell

)

Free Surface Processor 6: Mass Excess.

Definition at line 352 of file freeSurfacePostProcessor3D.hh.

                                                                          {
  using namespace olb::FreeSurface;
 
  // Note that EPSILON cannot be set in this operator because it is needed for the normal computation,
  // thus MASS is set here and EPSILON is set by the next operator.
  // Skip if the current cell is not an interface cell.
  if (!isCellType(cell, FreeSurface::Type::Interface)) {
    return;
  }
 
  T mass_excess = T(0);
  // If an interface cell is flagged as toGas,
  if (hasCellFlags(cell, FreeSurface::Flags::ToGas)) {
    // Get mass of the current interface cell
    const T mass = cell.template getField<FreeSurface::MASS>();
    mass_excess = mass;
    cell.template setField<FreeSurface::MASS>(T(0));
  }
  // If an interface cell is flagged as toFluid,
  else if (hasCellFlags(cell, FreeSurface::Flags::ToFluid)) {
    // Get density and mass of the current interface cell
    const T rho = cell.computeRho();
    const T mass = cell.template getField<FreeSurface::MASS>();
    mass_excess = mass - rho;
    cell.template setField<FreeSurface::MASS>(rho);
  }
  else { return; }
 
  // Redistribute the excess mass among the neighbouring interface cells, note that
  // if no neighbouring interface cells are found, the excess mass is lost or gained
  // by the current cell.
  // TODO: Thuerey Paper says we can't use new interface cells or flagged cells
  std::uint8_t oldIntefaceNbrs = std::uint8_t(0);
  std::uint8_t newIntefaceNbrs = std::uint8_t(0);
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto nbrCell = cell.neighbor(descriptors::c<DESCRIPTOR>(iPop));
 
    // Note that CELL_TYPE is not updated yet, and if no flag is set for this neighbour cell, it is an old interface cell
    // This means that this neighbour cell will remain an interface cell for the next time step
    if (isCellType(nbrCell, FreeSurface::Type::Interface) && hasCellFlags(nbrCell, FreeSurface::Flags::None)) {
      ++oldIntefaceNbrs;
    }
    // If the neighbour cell is flagged as Newinterface, it is a new interface cell
    else if (hasCellFlags(nbrCell, FreeSurface::Flags::NewInterface)) {
      ++newIntefaceNbrs;
    }
  }
 
  // Return if no interface cell is available to distribute the excess mass, i.e., the mass is lost or gained.
  const std::uint8_t intefaceNbrs = newIntefaceNbrs + oldIntefaceNbrs;
  if (intefaceNbrs == std::uint8_t(0)) { return; }
 
  // Distribution of excess mass among the neighbouring interface cells, weighted or evenly.
  bool enableAllInterfaces = true;
  constexpr bool weightedDistribution = true;
  Vector<T, DESCRIPTOR::q> weights{};
  T weightsSum = T(0);
 
  // (1) Distribute excess mass among the neighbouring interface cells, weighted by normal vector.
  if constexpr (weightedDistribution) {
    // Compute the weights for the mass excess distribution using normal vector of the current interface cell.
    auto normal = computeInterfaceNormal(cell);
    computeMassExcessWeights(cell, normal, enableAllInterfaces, weights);
 
    // Calculate sum of the computed weights
    for (const auto& weight : weights) { weightsSum += weight; }
 
    if (util::fabs(weightsSum) < zeroThreshold<T>()) {
      // (a) If no old interface cell is available in normal direction,
      //     fall back to weighted all interface cells distribution model.
      if (!enableAllInterfaces) {
        enableAllInterfaces = true;
        weightsSum = T(0);
        computeMassExcessWeights(cell, normal, enableAllInterfaces, weights);
        for (const auto& weight : weights) { weightsSum += weight; }
 
        // (b) If no interface cell is available in normal direction,
        //     fall back to evenly all interface cells distribution model.
        if (util::fabs(weightsSum) < zeroThreshold<T>()) {
          for (auto& weight : weights) { weight = T(1); }
          weightsSum = T(intefaceNbrs);
        }
      }
      // (c) If all interface cells are selected but none is available in normal direction,
      //     fall back to evenly all interface cells distribution model.
      else {
        for (auto& weight : weights) { weight = T(1); }
        weightsSum = T(intefaceNbrs);
      }
    }
  }
  // (2) Distribute excess mass among the neighbouring interface cells, evenly.
  else {
    // (a) Use old neighbouring interface cells to distribute excess mass.
    if (!enableAllInterfaces && oldIntefaceNbrs > std::uint8_t(0)) {
      for (auto& weight : weights) { weight = T(1); }
      weightsSum = T(oldIntefaceNbrs);
    }
    // (b) Use all neighbouring interface cells to distribute excess mass.
    else {
      enableAllInterfaces = true;
      for (auto& weight : weights) { weight = T(1); }
      weightsSum = T(intefaceNbrs);
    }
  }
 
  Vector<T, DESCRIPTOR::q> mass_exchange{};
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto nbrCell = cell.neighbor(descriptors::c<DESCRIPTOR>(iPop));
 
    if (isCellType(nbrCell, FreeSurface::Type::Interface) && hasCellFlags(nbrCell, FreeSurface::Flags::None)) {
      mass_exchange[iPop] = mass_excess * weights[iPop] / weightsSum;
    }
    else if (hasCellFlags(nbrCell, FreeSurface::Flags::NewInterface) && enableAllInterfaces) {
      mass_exchange[iPop] = mass_excess * weights[iPop] / weightsSum;
    }
    else {
      mass_exchange[iPop] = T(0);
    }
  }
 
  cell.template setField<FreeSurface::TEMP_MASS_EXCHANGE>(mass_exchange);
}

References olb::descriptors::c(), olb::util::fabs(), olb::FreeSurface::Interface, olb::FreeSurface::NewInterface, olb::FreeSurface::None, olb::FreeSurface::ToFluid, and olb::FreeSurface::ToGas.

Here is the call graph for this function:

getPriority()

template<typename T , typename DESCRIPTOR >

int olb::FreeSurfaceMassExcessPostProcessor3D< T, DESCRIPTOR >::getPriority ( ) const

inline

Definition at line 132 of file freeSurfacePostProcessor3D.h.

                          {
    return 5;
  }

Member Data Documentation

scope

template<typename T , typename DESCRIPTOR >

OperatorScope olb::FreeSurfaceMassExcessPostProcessor3D< T, DESCRIPTOR >::scope = OperatorScope::PerCell

staticconstexpr

Definition at line 130 of file freeSurfacePostProcessor3D.h.

---

The documentation for this class was generated from the following files:

• src/dynamics/freeSurfacePostProcessor3D.h
• src/dynamics/freeSurfacePostProcessor3D.hh