olb::FreeSurface Namespace Reference

Classes
struct	CELL_FLAGS
struct	CELL_TYPE
struct	DROP_ISOLATED_CELLS
struct	EPSILON
struct	HAS_SURFACE_TENSION
struct	LONELY_THRESHOLD
struct	MASS
struct	NeighbourInfo
struct	PREVIOUS_VELOCITY
struct	Stage0
struct	Stage1
struct	Stage2
struct	Stage3
struct	Stage4
struct	Stage5
struct	SURFACE_TENSION_PARAMETER
struct	TEMP_MASS_EXCHANGE
struct	TRANSITION

Enumerations
enum class	SolverType { rowPivotingLU , completePivotingLU , columnPivotingQR }
enum class	Type : std::uint8_t { Gas = 0 , Interface = 1 , Fluid = 2 , Solid = 4 }
enum class	Flags : std::uint8_t { None = 0 , ToGas = 1 , ToFluid = 2 , NewInterface = 4 }

Functions
template<typename T >
constexpr T	zeroThreshold ()
template<typename T >
bool	isNearZero (T a, T threshold)
template<typename T , typename DESCRIPTOR >
void	initialize (SuperLattice< T, DESCRIPTOR > &sLattice)
template<typename CELL >
bool	isCellType (CELL &cell, const FreeSurface::Type &type)
template<typename CELL >
void	setCellType (CELL &cell, const FreeSurface::Type &type)
template<typename CELL >
bool	hasNeighbour (CELL &cell, const FreeSurface::Type &type)
template<typename CELL >
bool	hasCellFlags (CELL &cell, const FreeSurface::Flags &flags)
template<typename CELL >
void	setCellFlags (CELL &cell, const FreeSurface::Flags &flags)
template<typename CELL >
bool	hasNeighbourFlags (CELL &cell, const FreeSurface::Flags &flags)
template<typename CELL >
NeighbourInfo	getNeighbourInfo (CELL &cell)
template<typename CELL , typename V >
bool	isHealthyInterface (CELL &cell)
template<typename CELL , typename V >
void	computeMassExcessWeights (CELL &cell, const Vector< V, CELL::descriptor_t::d > &normal, const bool &enableAllInterfaces, Vector< V, CELL::descriptor_t::q > &weights)
template<typename CELL , typename V >
V	getClampedEpsilon (const CELL &cell)
template<typename CELL , typename V >
V	getSmoothEpsilon (CELL &cell)
template<typename CELL , typename V >
V	getClampedSmoothEpsilon (const CELL &cell)
template<typename CELL , typename V >
Vector< V, CELL::descriptor_t::d >	computeInterfaceNormal (CELL &cell)
template<typename T , int N>
void	iterativeRefinement (const std::array< T, N *N > &A, const std::array< T, N *N > &M, std::array< T, N > &x, const std::array< T, N > &b, const int Nsol, int maxIter)
template<typename T , int N>
void	rowPivotingLUSolver (std::array< T, N *N > &M, std::array< T, N > &x, std::array< T, N > &b, const int Nsol)
template<typename T , int N>
void	completePivotingLUSolver (std::array< T, N *N > &M, std::array< T, N > &x, std::array< T, N > &b, const int Nsol)
template<typename T , int N>
void	columnPivotingQRSolver (std::array< T, N *N > &M, std::array< T, N > &x, std::array< T, N > &b, const int Nsol)
template<typename T >
T	plicCubeReduced (const T &volume, const Vector< T, 3 > &n)
template<typename T >
T	plicCube (const T &volume, const Vector< T, 3 > &n)
template<typename T >
T	plicCubeInverse (const T &d0, const Vector< T, 3 > &n)
template<typename CELL , typename V >
V	computeCurvaturePLIC2D (CELL &cell)
template<typename CELL , typename V >
V	computeCurvaturePLIC3D (CELL &cell)
template<typename CELL , typename V >
V	computeCurvatureFDM3D (CELL &cell)
template<typename CELL , typename V >
V	computeCurvature (CELL &cell)

Variables
constexpr SolverType	solver = SolverType::completePivotingLU

Enumeration Type Documentation

Flags

enum class olb::FreeSurface::Flags : std::uint8_t

strong

Enumerator
None
ToGas
ToFluid
NewInterface

Definition at line 72 of file freeSurfaceHelpers.h.

                 : std::uint8_t {
  None = 0,
  ToGas = 1,
  ToFluid = 2,
  NewInterface = 4
};

SolverType

enum class olb::FreeSurface::SolverType

strong

Enumerator
rowPivotingLU
completePivotingLU
columnPivotingQR

Definition at line 58 of file freeSurfaceHelpers.h.

                      {
  rowPivotingLU,
  completePivotingLU,
  columnPivotingQR
};

Type

enum class olb::FreeSurface::Type : std::uint8_t

strong

Enumerator
Gas
Interface
Fluid
Solid

Definition at line 65 of file freeSurfaceHelpers.h.

                : std::uint8_t {
  Gas = 0,
  Interface = 1,
  Fluid = 2,
  Solid = 4
};

Function Documentation

columnPivotingQRSolver()

template<typename T , int N>

void olb::FreeSurface::columnPivotingQRSolver	(	std::array< T, N *N > &	M,
		std::array< T, N > &	x,
		std::array< T, N > &	b,
		const int	Nsol )

Definition at line 448 of file freeSurfaceHelpers.hh.

                                                                                                         {
  // Permutation array to record column swaps.
  std::array<int, N> col_perm{};
  for (int i = 0; i < N; ++i) { col_perm[i] = i; }
 
  // Pre-computed column norm array, i.e., L2-normalization of each column.
  std::array<T, N> col_norm{}, col_norm_tmp{};
  for (int i = 0; i < Nsol; ++i) {
    T norm = T(0);
    for (int j = 0; j < Nsol; ++j) { norm += M[j * N + i] * M[j * N + i]; }
    col_norm[i] = util::sqrt(norm);
    col_norm_tmp[i] = col_norm[i];
  }
 
  // QR decomposition method with column pivoting.
  for (int i = 0; i < Nsol; ++i) {
    int pivotCol = i;
    T maxValue = col_norm[i];
 
    for (int j = i + 1; j < Nsol; ++j) {
      if (col_norm[j] > maxValue) {
        maxValue = col_norm[j];
        pivotCol = j;
      }
    }
 
    if (pivotCol != i) {
      for (int k = 0; k < Nsol; ++k) { std::swap(M[k * N + i], M[k * N + pivotCol]); }
      std::swap(col_perm[i], col_perm[pivotCol]);
      std::swap(col_norm[i], col_norm[pivotCol]);
      std::swap(col_norm_tmp[i], col_norm_tmp[pivotCol]);
    }
 
    // Compute the L2-norm of the i-th column, i.e., range is from row (i) to (Nsol - 1),
    // and the reflection factor (r)
    T r = T(0);
    {
      T norm = T(0);
      for (int k = i; k < Nsol; ++k) { norm += M[k * N + i] * M[k * N + i]; }
      norm = util::sqrt(norm);
      r = M[i * N + i] >= T(0) ? -norm : norm;
    }
 
    // Prepare and normalize the Householder array for the i-th column.
    std::array<T, N> v{};
    {
      v[0] = M[i * N + i] - r;
      for (int k = i + 1; k < Nsol; ++k) { v[k - i] = M[k * N + i]; }
 
      T norm = T(0);
      for (int k = 0; k < Nsol - i; ++k) { norm += v[k] * v[k]; }
      norm = util::sqrt(norm);
 
      // Otherwise, the Householder array is already zero and no reflection is applied.
      if (norm != T(0)) {
        for (int k = 0; k < Nsol - i; ++k) { v[k] /= norm; }
      }
    }
 
    // Apply the Householder reflection to the remaining columns.
    for (int k = i; k < Nsol; ++k) {
      T dot_product = T(0);
      for (int j = 0; j < Nsol - i; ++j) { dot_product += v[j] * M[(j + i) * N + k]; }
 
      dot_product *= T(2);
      for (int j = 0; j < Nsol - i; ++j) { M[(j + i) * N + k] -= v[j] * dot_product; }
    }
 
    // Apply the Householder reflection to the right-hand side array.
    {
      T dot_product = T(0);
      for (int k = 0; k < Nsol - i; ++k) { dot_product += v[k] * b[k + i]; }
 
      dot_product *= T(2);
      for (int k = 0; k < Nsol - i; ++k) { b[k + i] -= v[k] * dot_product; }
    }
 
    // Update the column norm array, i.e., adapted from LAPACK’s DGEQP3 routine.
    for (int k = i + 1; k < Nsol; ++k) {
      T tmp = col_norm[k] * col_norm[k] - M[i * N + k] * M[i * N + k];
 
      // Recompute L2-norm from scratch, if tmp becomes negative or it drops below
      // a specified threshold, e.g., T(0.1)
      if (tmp <= T(0) || util::sqrt(tmp) < T(0.1) * col_norm_tmp[k]) {
        T norm = T(0);
        for (int j = i + 1; j < Nsol; ++j) { norm += M[j * N + k] * M[j * N + k]; }
        col_norm[k] = util::sqrt(norm);
        col_norm_tmp[k] = col_norm[k];
      }
      else {
        col_norm[k] = util::sqrt(tmp);
      }
    }
 
    // For the i-th column, replace diagonal elements with (r) and the non-diagonal elements with zero.
    M[i * N + i] = r;
    for (int k = i + 1; k < Nsol; ++k) { M[k * N + i] = T(0); }
  }
 
  // Backward substitution step, i.e., solve R*x = Qᵀ*b in place.
  for (int i = Nsol - 1; i >= 0; --i) {
    T sum = T(0);
    for (int j = i + 1; j < Nsol; ++j) { sum += M[i * N + j] * x[j]; }
    x[i] = (b[i] - sum) / M[i * N + i];
  }
 
  // Reconstruct x by inverting the permutation, i.e., for each
  // index i, the original index is given by col_perm[i].
  std::array<T, N> tmp = x;
  for (int i = 0; i < Nsol; ++i) { x[col_perm[i]] = tmp[i]; }
}

References olb::norm(), and olb::util::sqrt().

Here is the call graph for this function:

completePivotingLUSolver()

template<typename T , int N>

void olb::FreeSurface::completePivotingLUSolver	(	std::array< T, N *N > &	M,
		std::array< T, N > &	x,
		std::array< T, N > &	b,
		const int	Nsol )

Definition at line 383 of file freeSurfaceHelpers.hh.

                                                                                                           {
  // Add Tikhonov Regularization to the diagonal elements of M
  for (int i = 0; i < Nsol; ++i) { M[i * N + i] += T(0); }
 
  // Permutation vectors to record row and columns swaps.
  std::array<int, N> row_perm{}, col_perm{};
  for (int i = 0; i < N; ++i) {
    row_perm[i] = i;
    col_perm[i] = i;
  }
 
  // LU decomposition method with complete pivoting.
  for (int i = 0; i < Nsol; ++i) {
    int pivotRow = i, pivotCol = i;
    T maxValue = std::abs(M[i * N + i]);
 
    for (int r = i; r < Nsol; ++r) {
      for (int c = i; c < Nsol; ++c) {
        T value = std::abs(M[r * N + c]);
        if (value > maxValue) {
          maxValue = value;
          pivotRow = r;
          pivotCol = c;
        }
      }
    }
 
    if (pivotRow != i) {
      for (int k = 0; k < Nsol; ++k) { std::swap(M[i * N + k], M[pivotRow * N + k]); }
      std::swap(row_perm[i], row_perm[pivotRow]);
      std::swap(b[i], b[pivotRow]);
    }
 
    if (pivotCol != i) {
      for (int k = 0; k < Nsol; ++k) { std::swap(M[k * N + i], M[k * N + pivotCol]); }
      std::swap(col_perm[i], col_perm[pivotCol]);
    }
 
    for (int j = i + 1; j < Nsol; ++j) {
      M[j * N + i] /= M[i * N + i];
      for (int k = i + 1; k < Nsol; ++k) { M[j * N + k] -= M[j * N + i] * M[i * N + k]; }
    }
  }
 
  // Copy the right-hand side array into x array
  for (int i = 0; i < Nsol; ++i) { x[i] = b[i]; }
 
  // Forward substitution step, i.e., solve L * x = b in place.
  for (int i = 0; i < Nsol; ++i) {
    for (int j = 0; j < i; ++j) { x[i] -= M[i * N + j] * x[j]; }
  }
 
  // Backward substitution step, i.e., solve U * x = b in place.
  for (int i = Nsol - 1; i >= 0; --i) {
    for (int j = i + 1; j < Nsol; ++j) { x[i] -= M[i * N + j] * x[j]; }
    x[i] /= M[i * N + i];
  }
 
  // Reconstruct original x by inverting the permutation, i.e.,
  // for each index i, the original index is given by col_perm[i].
  std::array<T, N> tmp = x;
  for (int i = 0; i < Nsol; ++i) { x[col_perm[i]] = tmp[i]; }
}

computeCurvature()

template<typename CELL , typename V >

V olb::FreeSurface::computeCurvature

(

CELL &

cell

)

Definition at line 910 of file freeSurfaceHelpers.hh.

                               {
  using DESCRIPTOR = typename CELL::descriptor_t;
 
  if constexpr (DESCRIPTOR::d == 2) {
    return computeCurvaturePLIC2D(cell);
  } else if (DESCRIPTOR::d == 3) {
    return computeCurvaturePLIC3D(cell);
  }
 
  return V(0);
}

computeCurvatureFDM3D()

template<typename CELL , typename V >

V olb::FreeSurface::computeCurvatureFDM3D

(

CELL &

cell

)

Definition at line 867 of file freeSurfaceHelpers.hh.

                                    {
  using DESCRIPTOR = typename CELL::descriptor_t;
 
  const auto normal = computeInterfaceNormal(cell);
  if (norm_squared(normal) < zeroThreshold<V>()) { return V(0); }
 
  const std::uint32_t gaussianMultipliers [DESCRIPTOR::q] =
  {
    std::uint32_t(8),
    std::uint32_t(4), std::uint32_t(4), std::uint32_t(4),
    std::uint32_t(2), std::uint32_t(2), std::uint32_t(2),
    std::uint32_t(2), std::uint32_t(2), std::uint32_t(2),
    std::uint32_t(1), std::uint32_t(1), std::uint32_t(1), std::uint32_t(1),
    std::uint32_t(4), std::uint32_t(4), std::uint32_t(4),
    std::uint32_t(2), std::uint32_t(2), std::uint32_t(2),
    std::uint32_t(2), std::uint32_t(2), std::uint32_t(2),
    std::uint32_t(1), std::uint32_t(1), std::uint32_t(1), std::uint32_t(1)
  };
 
  V curvature = V(0);
  V weightSum = V(0);
 
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto direction = descriptors::c<DESCRIPTOR>(iPop);
    auto nbrCell = cell.neighbor(direction);
 
    // Check if the neighbour is an interface cell or has a gas neighbour
    if (!isCellType(nbrCell, FreeSurface::Type::Interface) || !hasNeighbour(nbrCell, FreeSurface::Type::Gas)) {
      continue;
    }
 
    const Vector<V, DESCRIPTOR::d> ei = {V(direction[0]), V(direction[1]), V(direction[2])};
    const auto nbrNormal = util::normalize(computeInterfaceNormal(nbrCell));
    const V weight = gaussianMultipliers[iPop];
    weightSum += weight;
    curvature += weight * util::dotProduct(ei, nbrNormal);
  }
 
  curvature /= weightSum;
  return util::max(V(-1), util::min(V(1), curvature));
}

References olb::descriptors::c(), olb::util::dotProduct(), Gas, Interface, olb::util::max(), olb::util::min(), olb::norm_squared(), olb::util::normalize(), and zeroThreshold().

Here is the call graph for this function:

computeCurvaturePLIC2D()

template<typename CELL , typename V >

V olb::FreeSurface::computeCurvaturePLIC2D

(

CELL &

cell

)

Definition at line 644 of file freeSurfaceHelpers.hh.

                                     {
  using DESCRIPTOR = typename CELL::descriptor_t;
 
  // Setup a new coordinate system where bz is perpendicular to the surface, while bx and by are
  // tangent to the surface.
  const auto normal = computeInterfaceNormal(cell);
  Vector<V, 3> by = {normal[0], normal[1], V(0)};
  if (norm_squared(by) < zeroThreshold<V>()) { return V(0); }
  by = util::normalize(by);
 
  const Vector<V, 3> rn{V(0), V(0), V(1)};
  const Vector<V, 3> bx = crossProduct(by, rn);
 
  // Number of healthy interface neighbours
  std::uint8_t nbrInterfaces = std::uint8_t(0);
 
  // Compute z-offset of center point using PLIC cube
  const V centerOffset = plicCube<V>(getClampedEpsilon(cell), by);
 
  // Number of neighbouring interface points, i.e., less than or equal to 8 (neighbours) - 1 (liquid) - 1 (gas).
  std::array<Vector<V, DESCRIPTOR::d>, 6> points;
 
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto nbrCell = cell.neighbor(descriptors::c<DESCRIPTOR>(iPop));
 
    // Check if the neighbour is an interface cell or has a gas neighbour
    if (!isCellType(nbrCell, FreeSurface::Type::Interface) || !hasNeighbour(nbrCell, FreeSurface::Type::Gas)) {
      continue;
    }
 
    // Here it is assumed that neighbor normal vector is the same as center normal vector
    const auto direction = descriptors::c<DESCRIPTOR>(iPop);
    const Vector<V, 3> ei = {V(direction[0]), V(direction[1]), V(0)};
    const V offset = plicCube<V>(getClampedEpsilon(nbrCell), by) - centerOffset;
 
    // Transform coordinate system into (x, f(x)) and apply PLIC offset
    points[nbrInterfaces++] = {util::dotProduct(ei, bx), util::dotProduct(ei, by) + offset};
  }
 
  // Prepare the linear equation system, to be initialized with zeros
  std::array<V, 4> M{};
  std::array<V, 2> solution{};
  std::array<V, 2> b{};
 
  for (std::uint8_t i = 0; i < nbrInterfaces; ++i) {
    const V x = points[i][0];
    const V y = points[i][1];
    const V x2 = util::sqr(x);
    const V x3 = util::sqr(x) * x;
 
    M[0] += x2 * x2;
    M[1] += x3;
    M[3] += x2;
 
    b[0] += x2 * y;
    b[1] += x * y;
  }
 
  // Fill the lower triangle of the symmetric matrix
  M[2] = M[1];
 
  // Solve the linear system of equations, compile-time selection.
  if constexpr (solver == SolverType::rowPivotingLU) {
    if (nbrInterfaces >= 2) {
      rowPivotingLUSolver<V, 2>(M, solution, b, 2);
    } else {
      rowPivotingLUSolver<V, 2>(M, solution, b, util::min(2, nbrInterfaces));
    }
  }
  else if constexpr (solver == SolverType::completePivotingLU) {
    if (nbrInterfaces >= 2) {
      completePivotingLUSolver<V, 2>(M, solution, b, 2);
    } else {
      completePivotingLUSolver<V, 2>(M, solution, b, util::min(2, nbrInterfaces));
    }
  }
  else if constexpr (solver == SolverType::columnPivotingQR) {
    if (nbrInterfaces >= 2) {
      columnPivotingQRSolver<V, 2>(M, solution, b, 2);
    } else {
      columnPivotingQRSolver<V, 2>(M, solution, b, util::min(2, nbrInterfaces));
    }
  }
  else {
    // Fall back to default solver
    if (nbrInterfaces >= 2) {
      completePivotingLUSolver<V, 2>(M, solution, b, 2);
    } else {
      completePivotingLUSolver<V, 2>(M, solution, b, util::min(2, nbrInterfaces));
    }
  }
 
  const V A = solution[0], H = solution[1];
  const V curvature = V(2) * A * util::cube(V(1) / util::sqrt(H * H + V(1)));
  return util::max(V(-1), util::min(V(1), curvature));
}

References olb::descriptors::c(), columnPivotingQR, completePivotingLU, olb::crossProduct(), olb::util::cube(), olb::util::dotProduct(), Gas, Interface, olb::util::max(), olb::util::min(), olb::norm_squared(), olb::util::normalize(), rowPivotingLU, solver, olb::util::sqr(), olb::util::sqrt(), and zeroThreshold().

Here is the call graph for this function:

computeCurvaturePLIC3D()

template<typename CELL , typename V >

V olb::FreeSurface::computeCurvaturePLIC3D

(

CELL &

cell

)

Definition at line 743 of file freeSurfaceHelpers.hh.

                                     {
  using DESCRIPTOR = typename CELL::descriptor_t;
 
  // Setup a new coordinate system where bz is perpendicular to the surface, while bx and by are
  // tangent to the surface.
  auto bz = computeInterfaceNormal(cell);
  if (norm_squared(bz) < zeroThreshold<V>()) { return V(0); }
  bz = util::normalize(bz);
 
  // A random normalized vector that is just by random chance not collinear with bz
  // const Vector<V, DESCRIPTOR::d> rn{V(0.56270900), V(0.32704452), V(0.75921047)};
  const Vector<V, DESCRIPTOR::d> rn{V(0.5402151198529208), V(0.2765640977634561), V(0.7947829415070379)};
 
  // Normalization is necessary here because bz and rn are not perpendicular
  const Vector<V, DESCRIPTOR::d> by = util::normalize(crossProduct(bz, rn));
  const Vector<V, DESCRIPTOR::d> bx = crossProduct(by, bz);
 
  // Number of healthy interface neighbours
  std::uint8_t nbrInterfaces = std::uint8_t(0);
 
  // Compute z-offset of center point using PLIC cube
  const V centerOffset = plicCube<V>(getClampedEpsilon(cell), bz);
 
  // Number of neighbouring interface points, i.e., less than or equal to 26 (neighbours) - 1 (liquid) - 1 (gas).
  std::array<Vector<V, DESCRIPTOR::d>, 24> points;
 
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto nbrCell = cell.neighbor(descriptors::c<DESCRIPTOR>(iPop));
 
    // Check if the neighbour is an interface cell or has a gas neighbour
    if (!isCellType(nbrCell, FreeSurface::Type::Interface) || !hasNeighbour(nbrCell, FreeSurface::Type::Gas)) {
      continue;
    }
 
    // Here it is assumed that neighbor normal vector is the same as center normal vector
    const auto direction = descriptors::c<DESCRIPTOR>(iPop);
    const Vector<V, DESCRIPTOR::d> ei = {V(direction[0]), V(direction[1]), V(direction[2])};
    const V offset = plicCube<V>(getClampedEpsilon(nbrCell), bz) - centerOffset;
 
    // Transform coordinate system into (x, y, f(x,y)) and apply PLIC offset
    points[nbrInterfaces++] = {util::dotProduct(ei, bx), util::dotProduct(ei, by), util::dotProduct(ei, bz) + offset};
  }
 
  // Prepare the linear equation system, to be initialized with zeros
  std::array<V, 25> M{};
  std::array<V, 5> solution{};
  std::array<V, 5> b{};
 
  for (std::uint8_t i = 0; i < nbrInterfaces; ++i) {
    const V x = points[i][0];
    const V y = points[i][1];
    const V z = points[i][2];
    const V x2 = util::sqr(x);
    const V y2 = util::sqr(y);
    const V x3 = x2 * x;
    const V y3 = y2 * y;
 
    M[0] += x2 * x2;
    M[1] += x2 * y2;
    M[2] += x3 * y;
    M[3] += x3;
    M[4] += x2 * y;
    M[6] += y2 * y2;
    M[7] += x * y3;
    M[8] += x * y2;
    M[9] += y3;
    M[12] += x2 * y2;
    M[13] += x2 * y;
    M[14] += x * y2;
    M[18] += x2;
    M[19] += x * y;
    M[24] += y2;
 
    b[0] += x2 * z;
    b[1] += y2 * z;
    b[2] += x * y * z;
    b[3] += x * z;
    b[4] += y * z;
  }
 
  // Fill the lower triangle of the symmetric matrix
  for (std::uint8_t i = 1; i < 5; ++i) {
    for (std::uint8_t j = 0; j < i; ++j) {
      M[i * 5 + j] = M[j * 5 + i];
    }
  }
 
  // Solve the linear system of equations, compile-time selection.
  if constexpr (solver == SolverType::rowPivotingLU) {
    if (nbrInterfaces >= 5) {
      rowPivotingLUSolver<V, 5>(M, solution, b, 5);
    } else {
      rowPivotingLUSolver<V, 5>(M, solution, b, util::min(5, nbrInterfaces));
    }
  }
  else if constexpr (solver == SolverType::completePivotingLU) {
    if (nbrInterfaces >= 5) {
      completePivotingLUSolver<V, 5>(M, solution, b, 5);
    } else {
      completePivotingLUSolver<V, 5>(M, solution, b, util::min(5, nbrInterfaces));
    }
  }
  else if constexpr (solver == SolverType::columnPivotingQR) {
    if (nbrInterfaces >= 5) {
      columnPivotingQRSolver<V, 5>(M, solution, b, 5);
    } else {
      columnPivotingQRSolver<V, 5>(M, solution, b, util::min(5, nbrInterfaces));
    }
  }
  else {
    // Fall back to default solver
    if (nbrInterfaces >= 5) {
      completePivotingLUSolver<V, 5>(M, solution, b, 5);
    } else {
      completePivotingLUSolver<V, 5>(M, solution, b, util::min(5, nbrInterfaces));
    }
  }
 
  const V A = solution[0], B = solution[1], C = solution[2], H = solution[3], I = solution[4];
  const V curvature = (A * (I * I + V(1)) + B * (H * H + V(1)) - C * H * I) * util::cube(V(1) / util::sqrt(H * H + I * I + V(1)));
  return util::max(V(-1), util::min(V(1), curvature));
}

References olb::descriptors::c(), columnPivotingQR, completePivotingLU, olb::crossProduct(), olb::util::cube(), olb::util::dotProduct(), Gas, Interface, olb::util::max(), olb::util::min(), olb::norm_squared(), olb::util::normalize(), rowPivotingLU, solver, olb::util::sqr(), olb::util::sqrt(), and zeroThreshold().

Here is the call graph for this function:

computeInterfaceNormal()

template<typename CELL , typename V >

Vector< V, CELL::descriptor_t::d > olb::FreeSurface::computeInterfaceNormal

(

CELL &

cell

)

Definition at line 223 of file freeSurfaceHelpers.hh.

                                                                {
 
  // Compute interface normal using Parker-Young approximation
  Vector<V, CELL::descriptor_t::d> normal{};
  if constexpr (CELL::descriptor_t::d == 3) {
    auto eps_minusY             = getClampedEpsilon(cell.neighbor({0, -1, 0}));
    auto eps_plusY              = getClampedEpsilon(cell.neighbor({0, 1, 0}));
    auto eps_minusX             = getClampedEpsilon(cell.neighbor({-1, 0, 0}));
    auto eps_plusX              = getClampedEpsilon(cell.neighbor({1, 0, 0}));
    auto eps_minusYminusX       = getClampedEpsilon(cell.neighbor({-1, -1, 0}));
    auto eps_minusYplusX        = getClampedEpsilon(cell.neighbor({1, -1, 0}));
    auto eps_plusYminusX        = getClampedEpsilon(cell.neighbor({-1, 1, 0}));
    auto eps_plusYplusX         = getClampedEpsilon(cell.neighbor({1, 1, 0}));
    auto eps_minusZ             = getClampedEpsilon(cell.neighbor({0, 0, -1}));
    auto eps_plusZ              = getClampedEpsilon(cell.neighbor({0, 0, 1}));
    auto eps_minusZminusY       = getClampedEpsilon(cell.neighbor({0, -1, -1}));
    auto eps_minusZplusY        = getClampedEpsilon(cell.neighbor({0, 1, -1}));
    auto eps_minusZminusX       = getClampedEpsilon(cell.neighbor({-1, 0, -1}));
    auto eps_minusZplusX        = getClampedEpsilon(cell.neighbor({1, 0, -1}));
    auto eps_plusZminusY        = getClampedEpsilon(cell.neighbor({0, -1, 1}));
    auto eps_plusZplusY         = getClampedEpsilon(cell.neighbor({0, 1, 1}));
    auto eps_plusZminusX        = getClampedEpsilon(cell.neighbor({-1, 0, 1}));
    auto eps_plusZplusX         = getClampedEpsilon(cell.neighbor({1, 0, 1}));
    auto eps_minusZminusYminusX = getClampedEpsilon(cell.neighbor({-1, -1, -1}));
    auto eps_minusZplusYminusX  = getClampedEpsilon(cell.neighbor({-1, 1, -1}));
    auto eps_minusZminusYplusX  = getClampedEpsilon(cell.neighbor({1, -1, -1}));
    auto eps_minusZplusYplusX   = getClampedEpsilon(cell.neighbor({1, 1, -1}));
    auto eps_plusZminusYminusX  = getClampedEpsilon(cell.neighbor({-1, -1, 1}));
    auto eps_plusZplusYminusX   = getClampedEpsilon(cell.neighbor({-1, 1, 1}));
    auto eps_plusZminusYplusX   = getClampedEpsilon(cell.neighbor({1, -1, 1}));
    auto eps_plusZplusYplusX    = getClampedEpsilon(cell.neighbor({1, 1, 1}));
 
    // Rearrange arithmetic operations to minimize the number of cancellation events, i.e., positive and negative terms.
    normal[0]  = V(4) * (eps_minusX - eps_plusX);
    normal[0] += V(2) * ((eps_minusYminusX + eps_minusZminusX + eps_plusYminusX + eps_plusZminusX) - (eps_plusYplusX + eps_plusZplusX + eps_minusYplusX + eps_minusZplusX));
    normal[0] += V(1) * ((eps_minusZminusYminusX + eps_plusZminusYminusX + eps_minusZplusYminusX + eps_plusZplusYminusX) - (eps_plusZplusYplusX + eps_minusZplusYplusX + eps_plusZminusYplusX + eps_minusZminusYplusX));
 
    normal[1]  = V(4) * (eps_minusY - eps_plusY);
    normal[1] += V(2) * ((eps_minusYminusX + eps_minusZminusY + eps_minusYplusX + eps_plusZminusY) - (eps_plusYplusX + eps_plusZplusY + eps_plusYminusX + eps_minusZplusY));
    normal[1] += V(1) * ((eps_minusZminusYminusX + eps_plusZminusYminusX + eps_plusZminusYplusX + eps_minusZminusYplusX) - (eps_plusZplusYplusX + eps_minusZplusYplusX + eps_minusZplusYminusX + eps_plusZplusYminusX));
 
    normal[2]  = V(4) * (eps_minusZ - eps_plusZ);
    normal[2] += V(2) * ((eps_minusZminusX + eps_minusZminusY + eps_minusZplusX + eps_minusZplusY) - (eps_plusZplusX + eps_plusZplusY + eps_plusZminusX + eps_plusZminusY));
    normal[2] += V(1) * ((eps_minusZminusYminusX + eps_minusZplusYplusX + eps_minusZplusYminusX + eps_minusZminusYplusX) - (eps_plusZplusYplusX + eps_plusZminusYminusX + eps_plusZminusYplusX + eps_plusZplusYminusX));
 
    return normal;
  }
  else {
    auto eps_minusY       = getClampedEpsilon(cell.neighbor({0, -1}));
    auto eps_plusY        = getClampedEpsilon(cell.neighbor({0, 1}));
    auto eps_minusX       = getClampedEpsilon(cell.neighbor({-1, 0}));
    auto eps_plusX        = getClampedEpsilon(cell.neighbor({1, 0}));
    auto eps_minusYminusX = getClampedEpsilon(cell.neighbor({-1, -1}));
    auto eps_minusYplusX  = getClampedEpsilon(cell.neighbor({1, -1}));
    auto eps_plusYminusX  = getClampedEpsilon(cell.neighbor({-1, 1}));
    auto eps_plusYplusX   = getClampedEpsilon(cell.neighbor({1, 1}));
 
    // Rearrange arithmetic operations to minimize the number of cancellation events, i.e., positive and negative terms.
    normal[0] = V(2) * (eps_minusX - eps_plusX);
    normal[0] += V(1) * ((eps_minusYminusX + eps_plusYminusX) - (eps_plusYplusX + eps_minusYplusX));
 
    normal[1] = V(2) * (eps_minusY - eps_plusY);
    normal[1] += V(1) * ((eps_minusYminusX + eps_minusYplusX) - (eps_plusYplusX + eps_plusYminusX));
 
    return normal;
  }
}

computeMassExcessWeights()

template<typename CELL , typename V >

void olb::FreeSurface::computeMassExcessWeights	(	CELL &	cell,
		const Vector< V, CELL::descriptor_t::d > &	normal,
		const bool &	enableAllInterfaces,
		Vector< V, CELL::descriptor_t::q > &	weights )

Definition at line 131 of file freeSurfaceHelpers.hh.

  {
  using DESCRIPTOR = typename CELL::descriptor_t;
 
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto direction = descriptors::c<DESCRIPTOR>(iPop);
    auto nbrCell = cell.neighbor(direction);
 
    Vector<V, DESCRIPTOR::d> ei{};
    if constexpr (DESCRIPTOR::d == 3) {
      ei = {V(direction[0]), V(direction[1]), V(direction[2])};
    }
    else { ei = {V(direction[0]), V(direction[1])}; }
 
    if (isCellType(nbrCell, FreeSurface::Type::Interface) && hasCellFlags(nbrCell, FreeSurface::Flags::None)) {
      const V nDotDirection = util::dotProduct(normal, ei);
 
      if (hasCellFlags(cell, FreeSurface::Flags::ToFluid)) {
        // This cell was converted from interface to fluid, so the normal vector is used as is
        weights[iPop] = nDotDirection > V(0) ? nDotDirection : V(0);
      }
      else if (hasCellFlags(cell, FreeSurface::Flags::ToGas)) {
        // This cell was converted from interface to gas, so the normal vector is inverted
        weights[iPop] = nDotDirection < V(0) ? -nDotDirection : V(0);
      }
    }
    else if (hasCellFlags(nbrCell, FreeSurface::Flags::NewInterface) && enableAllInterfaces) {
      const V nDotDirection = util::dotProduct(normal, ei);
 
      if (hasCellFlags(cell, FreeSurface::Flags::ToFluid)) {
        // This cell was converted from interface to fluid, so the normal vector is used as is
        weights[iPop] = nDotDirection > V(0) ? nDotDirection : V(0);
      }
      else if (hasCellFlags(cell, FreeSurface::Flags::ToGas)) {
        // This cell was converted from interface to gas, so the normal vector is inverted
        weights[iPop] = nDotDirection < V(0) ? -nDotDirection : V(0);
      }
    }
    else {
      // If the neighbour is not an interface cell, the weight is set to zero
      weights[iPop] = V(0);
    }
  }
}

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

Here is the call graph for this function:

getClampedEpsilon()

template<typename CELL , typename V >

V olb::FreeSurface::getClampedEpsilon

(

const CELL &

cell

)

Definition at line 180 of file freeSurfaceHelpers.hh.

                                      {
 
  V epsilon = cell.template getField<FreeSurface::EPSILON>();
  return util::max(V(0), util::min(V(1), epsilon));
}

References olb::util::max(), and olb::util::min().

Here is the call graph for this function:

getClampedSmoothEpsilon()

template<typename CELL , typename V >

V olb::FreeSurface::getClampedSmoothEpsilon

(

const CELL &

cell

)

Definition at line 216 of file freeSurfaceHelpers.hh.

                                            {
 
  V epsilon = getSmoothEpsilon(cell);
  return util::max(V(0), util::min(V(1), epsilon));
}

References olb::util::max(), and olb::util::min().

Here is the call graph for this function:

getNeighbourInfo()

template<typename CELL >

NeighbourInfo olb::FreeSurface::getNeighbourInfo

(

CELL &

cell

)

Definition at line 86 of file freeSurfaceHelpers.hh.

                                           {
  NeighbourInfo info{};
  using DESCRIPTOR = typename CELL::descriptor_t;
 
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto nbrCell = cell.neighbor(descriptors::c<DESCRIPTOR>(iPop));
 
    if (isCellType(nbrCell, FreeSurface::Type::Gas)) {
      info.has_gas_neighbours = true;
    }
    else if (isCellType(nbrCell, FreeSurface::Type::Fluid)) {
      info.has_fluid_neighbours = true;
    }
    else if (isCellType(nbrCell, FreeSurface::Type::Interface)) {
      ++info.interface_neighbours;
    }
  }
 
  return info;
}

References olb::descriptors::c(), Fluid, Gas, and Interface.

Here is the call graph for this function:

getSmoothEpsilon()

template<typename CELL , typename V >

V olb::FreeSurface::getSmoothEpsilon

(

CELL &

cell

)

Definition at line 187 of file freeSurfaceHelpers.hh.

                               {
  using DESCRIPTOR = typename CELL::descriptor_t;
 
  // K8 kernel of Williams et al.
  // Avoid uisng expensive std::pow to boost performance
  auto kernel = [](auto radius, const auto& direction) {
    auto norm = norm_squared(direction);
    if (norm < radius * radius) {
      auto tmp = decltype(radius)(1) - norm / (radius * radius);
      return tmp * tmp * tmp * tmp;
    } else { return decltype(radius)(0); }
  };
 
  // Kernel support radius is fixed at 2.0
  V radius = V(2);
  V denominator = V(0);
  V tmp = V(0);
 
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto nbrCell = cell.neighbor(descriptors::c<DESCRIPTOR>(iPop));
    const auto direction = descriptors::c<DESCRIPTOR>(iPop);
    tmp += kernel(radius, direction) * nbrCell.template getField<FreeSurface::EPSILON>();
    if (norm(direction) < radius) { denominator += kernel(radius, direction); }
  }
 
  return tmp / denominator;
}

References olb::descriptors::c(), olb::norm(), and olb::norm_squared().

Here is the call graph for this function:

hasCellFlags()

template<typename CELL >

bool olb::FreeSurface::hasCellFlags	(	CELL &	cell,
		const FreeSurface::Flags &	flags )

Definition at line 64 of file freeSurfaceHelpers.hh.

                                                             {
  return static_cast<bool>(cell.template getField<FreeSurface::CELL_FLAGS>() & flags);
}

hasNeighbour()

template<typename CELL >

bool olb::FreeSurface::hasNeighbour	(	CELL &	cell,
		const FreeSurface::Type &	type )

Definition at line 52 of file freeSurfaceHelpers.hh.

                                                           {
  using DESCRIPTOR = typename CELL::descriptor_t;
 
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto nbrCell = cell.neighbor(descriptors::c<DESCRIPTOR>(iPop));
    if (isCellType(nbrCell, type)) { return true; }
  }
 
  return false;
}

References olb::descriptors::c().

Here is the call graph for this function:

hasNeighbourFlags()

template<typename CELL >

bool olb::FreeSurface::hasNeighbourFlags	(	CELL &	cell,
		const FreeSurface::Flags &	flags )

Definition at line 74 of file freeSurfaceHelpers.hh.

                                                                  {
  using DESCRIPTOR = typename CELL::descriptor_t;
 
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto nbrCell = cell.neighbor(descriptors::c<DESCRIPTOR>(iPop));
    if (hasCellFlags(nbrCell, flags)) { return true; }
  }
 
  return false;
}

References olb::descriptors::c().

Here is the call graph for this function:

initialize()

template<typename T , typename DESCRIPTOR >

void olb::FreeSurface::initialize

(

SuperLattice< T, DESCRIPTOR > &

sLattice

)

Definition at line 32 of file freeSurfaceHelpers.hh.

                                                      {
  sLattice.executePostProcessors(FreeSurface::Stage0());
  sLattice.executePostProcessors(FreeSurface::Stage1());
  sLattice.executePostProcessors(FreeSurface::Stage2());
  sLattice.executePostProcessors(FreeSurface::Stage3());
  sLattice.executePostProcessors(FreeSurface::Stage4());
  sLattice.executePostProcessors(FreeSurface::Stage5());
}

References olb::SuperLattice< T, DESCRIPTOR >::executePostProcessors().

Here is the call graph for this function:

isCellType()

template<typename CELL >

bool olb::FreeSurface::isCellType	(	CELL &	cell,
		const FreeSurface::Type &	type )

Definition at line 42 of file freeSurfaceHelpers.hh.

                                                         {
  return cell.template getField<FreeSurface::CELL_TYPE>() == type;
}

isHealthyInterface()

template<typename CELL , typename V >

bool olb::FreeSurface::isHealthyInterface

(

CELL &

cell

)

Definition at line 108 of file freeSurfaceHelpers.hh.

                                    {
  using DESCRIPTOR = typename CELL::descriptor_t;
  bool has_fluid_neighbours = false;
  bool has_gas_neighbours = false;
 
  if (!isCellType(cell, FreeSurface::Type::Interface)) { return false; }
 
  for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
    auto nbrCell = cell.neighbor(descriptors::c<DESCRIPTOR>(iPop));
    if (isCellType(nbrCell, FreeSurface::Type::Gas)) {
      has_gas_neighbours = true;
      if (has_fluid_neighbours) { return true; }
    }
    else if (isCellType(nbrCell, FreeSurface::Type::Fluid)) {
      has_fluid_neighbours = true;
      if (has_gas_neighbours) { return true; }
    }
  }
 
  return false;
}

References olb::descriptors::c(), Fluid, Gas, and Interface.

Here is the call graph for this function:

isNearZero()

template<typename T >

bool olb::FreeSurface::isNearZero ( T a, T threshold )

inline

Definition at line 54 of file freeSurfaceHelpers.h.

                                         {
  return util::abs(a) < threshold;
}

References olb::util::abs().

Here is the call graph for this function:

iterativeRefinement()

template<typename T , int N>

void olb::FreeSurface::iterativeRefinement	(	const std::array< T, N *N > &	A,
		const std::array< T, N *N > &	M,
		std::array< T, N > &	x,
		const std::array< T, N > &	b,
		const int	Nsol,
		int	maxIter )

Definition at line 292 of file freeSurfaceHelpers.hh.

  {
  std::array<T, N> residual{};
  for (int iter = 0; iter < maxIter; ++iter) {
    // Compute the residual, i.e., r = b - A * x
    for (int i = 0; i < Nsol; ++i) {
      T sum = T(0);
      for (int j = 0; j < Nsol; ++j) { sum += A[i * N + j] * x[j]; }
      residual[i] = b[i] - sum;
    }
 
    // Forward substitution
    std::array<T, N> dx = residual;
    for (int i = 0; i < Nsol; ++i) {
      for (int k = 0; k < i; ++k) { dx[i] -= M[i * N + k] * dx[k]; }
    }
 
    // Backward substitution
    for (int i = Nsol - 1; i >= 0; --i) {
      for (int k = i + 1; k < Nsol; ++k) { dx[i] -= M[i * N + k] * dx[k]; }
      dx[i] /= M[i * N + i];
    }
 
    // Update the solution, i.e., x = x + dx
    for (int i = 0; i < Nsol; ++i) { x[i] += dx[i]; }
  }
}

plicCube()

template<typename T >

T olb::FreeSurface::plicCube	(	const T &	volume,
		const Vector< T, 3 > &	n )

Definition at line 596 of file freeSurfaceHelpers.hh.

                                                   {
  // Unit cube-plane intersection: using volume and normal vector to compute the offset
  const T ax = util::abs(n[0]), ay = util::abs(n[1]), az = util::abs(n[2]);
 
  // Eliminate symmetry cases, and normalize n using L1 norm
  const T V = T(0.5) - util::abs(volume - T(0.5)), l = ax + ay + az;
  const T n1 = util::min(util::min(ax, ay), az) / l;
  const T n3 = util::max(util::max(ax, ay), az) / l;
  const T n2 = std::fdim(T(1), n1 + n3); // ensure n2 >= 0
  const T d = plicCubeReduced<T>(V, {n1, n2, n3});
 
  return l * std::copysign(T(0.5) - d, volume - T(0.5));
}

References olb::util::abs(), olb::util::max(), and olb::util::min().

Here is the call graph for this function:

plicCubeInverse()

template<typename T >

T olb::FreeSurface::plicCubeInverse	(	const T &	d0,
		const Vector< T, 3 > &	n )

Definition at line 612 of file freeSurfaceHelpers.hh.

                                                      {
  // Eliminate majority of cases due to symmetry
  const T n1 = util::min(util::min(util::abs(n[0]), util::abs(n[1])), util::abs(n[2]));
  const T n3 = util::max(util::max(util::abs(n[0]), util::abs(n[1])), util::abs(n[2]));
  const T n2 = util::abs(n[0]) - n1 + util::abs(n[1]) + util::abs(n[2]) - n3;
 
  // Compute PLIC using reduced symmetry, shifting origin from (0.0, 0.0, 0.0) to (0.5, 0.5, 0.5)
  const T d = T(0.5) * (n1 + n2 + n3) - util::abs(d0);
 
  T V = T(0);
  if (util::min(n1 + n2, n3) <= d && d <= n3) {
    // (I) Case 5
    V = (d - T(0.5) * (n1 + n2)) / n3;
  }
  else if (d < n1) {
    // (II) Case 1
    V = util::cube(d) / (T(6) * n1 * n2 * n3);
  }
  else if (d <= n2) {
    // (III) Case 2
    V = (T(3) * d * (d - n1) + util::sqr(n1)) / (T(6) * n2 * n3);
  }
  else {
    // (IV) Case 3 or 4
    V = (util::cube(d) - util::cube(d - n1) - util::cube(d - n2) - util::cube(std::fdim(d, n3))) / (T(6) * n1 * n2 * n3);
  }
 
  return std::copysign(T(0.5) - V, d0) + V(0.5);
}

References olb::util::abs(), olb::util::cube(), olb::util::max(), olb::util::min(), and olb::util::sqr().

Here is the call graph for this function:

plicCubeReduced()

template<typename T >

T olb::FreeSurface::plicCubeReduced	(	const T &	volume,
		const Vector< T, 3 > &	n )

Definition at line 562 of file freeSurfaceHelpers.hh.

                                                          {
  const T n1 = n[0], n2 = n[1], n3 = n[2];
  const T n12 = n1 + n2, n3V = n3 * volume;
 
  // (I) Case 5
  if (n12 <= T(2) * n3V) { return n3V + T(0.5) * n12; }
 
  // (II) After case 5, check n2 > 0 is true
  const T sqn1 = util::sqr(n1), n26 = T(6) * n2, v1 = sqn1 / n26;
 
  // (III) Case 2
  if (v1 <= n3V && n3V < v1 + T(0.5) * (n2 - n1)) {
    return T(0.5) * (n1 + util::sqrt(sqn1 + T(8) * n2 * (n3V - v1)));
  }
 
  // (IV) Case 1
  const T V6 = n1 * n26 * n3V;
  if (n3V < v1) { return std::cbrt(V6); }
 
  // (V) After case 2, check n1 > 0 is true
  const T v3 = n3 < n12 ? (util::sqr(n3) * (T(3) * n12 - n3) + sqn1 * (n1 - T(3) * n3) + util::sqr(n2) * (n2 - T(3) * n3)) / (n1 * n26) : T(0.5) * n12;
  const T sqn12 = sqn1 + util::sqr(n2), V6cbn12 = V6 - util::cube(n1) - util::cube(n2);
  const bool case34 = n3V < v3; // true: case (3), false: case (4)
  const T a = case34 ? V6cbn12 : T(0.5) * (V6cbn12 - util::cube(n3));
  const T b = case34 ? sqn12 : T(0.5) * (sqn12 + util::sqr(n3));
  const T c = case34 ? n12 : T(0.5);
  const T t = util::sqrt(util::sqr(c) - b);
  static constexpr T oneThird = T(1) / T(3);
 
  return c - T(2) * t * util::sin(oneThird * util::asin((util::cube(c) - T(0.5) * a - T(1.5) * b * c) / util::cube(t)));
}

References olb::util::asin(), olb::util::cube(), olb::util::sin(), olb::util::sqr(), and olb::util::sqrt().

Here is the call graph for this function:

rowPivotingLUSolver()

template<typename T , int N>

void olb::FreeSurface::rowPivotingLUSolver	(	std::array< T, N *N > &	M,
		std::array< T, N > &	x,
		std::array< T, N > &	b,
		const int	Nsol )

Definition at line 324 of file freeSurfaceHelpers.hh.

                                                                                                      {
  // Add Tikhonov Regularization to the diagonal elements of M
  for (int i = 0; i < Nsol; ++i) { M[i * N + i] += T(0); }
 
  // Make a copy of the original coefficient matrix.
  std::array<T, N * N> Prev_M = M;
 
  // Permutation vector to record row swaps.
  std::array<int, N> row_perm{};
  for (int i = 0; i < N; ++i) { row_perm[i] = i; }
 
  // LU decomposition method with row pivoting.
  for (int i = 0; i < Nsol; ++i) {
    int pivotRow = i;
    T maxValue = std::abs(M[i * N + i]);
 
    for (int j = i + 1; j < Nsol; ++j) {
      T value = std::abs(M[j * N + i]);
      if (value > maxValue) {
        maxValue = value;
        pivotRow = j;
      }
    }
 
    if (pivotRow != i) {
      for (int k = 0; k < Nsol; ++k) { std::swap(M[i * N + k], M[pivotRow * N + k]); }
      std::swap(row_perm[i], row_perm[pivotRow]);
      std::swap(b[i], b[pivotRow]);
    }
 
    for (int j = i + 1; j < Nsol; ++j) {
      M[j * N + i] /= M[i * N + i];
      for (int k = i + 1; k < Nsol; ++k) { M[j * N + k] -= M[j * N + i] * M[i * N + k]; }
    }
  }
 
  // Copy the right-hand side array into x array
  for (int i = 0; i < Nsol; ++i) { x[i] = b[i]; }
 
  // Forward substitution step, i.e., solve L * x = b in place.
  for (int i = 0; i < Nsol; ++i) {
    for (int j = 0; j < i; ++j) { x[i] -= M[i * N + j] * x[j]; }
  }
 
  // Backward substitution step, i.e., solve U * x = b in place.
  for (int i = Nsol - 1; i >= 0; --i) {
    for (int j = i + 1; j < Nsol; ++j) { x[i] -= M[i * N + j] * x[j]; }
    x[i] /= M[i * N + i];
  }
 
  // Perform iterative refinement of the solution.
  std::array<T, N * N> tmp = Prev_M;
  for (int i = 0; i < Nsol; ++i) {
    for (int j = 0; j < Nsol; ++j) { Prev_M[i * N + j] = tmp[row_perm[i] * N + j]; }
  }
  iterativeRefinement<T, N>(Prev_M, M, x, b, Nsol);
}

setCellFlags()

template<typename CELL >

void olb::FreeSurface::setCellFlags	(	CELL &	cell,
		const FreeSurface::Flags &	flags )

Definition at line 69 of file freeSurfaceHelpers.hh.

                                                             {
  cell.template setField<FreeSurface::CELL_FLAGS>(flags);
}

setCellType()

template<typename CELL >

void olb::FreeSurface::setCellType	(	CELL &	cell,
		const FreeSurface::Type &	type )

Definition at line 47 of file freeSurfaceHelpers.hh.

                                                          {
  cell.template setField<FreeSurface::CELL_TYPE>(type);
}

zeroThreshold()

template<typename T >

T olb::FreeSurface::zeroThreshold ( )

constexpr

Definition at line 45 of file freeSurfaceHelpers.h.

                            {
  if constexpr (std::is_same<T, float>::value) {
    return T{1e-6};
  } else {
    return T{1e-14};
  }
}

Here is the caller graph for this function:

Variable Documentation

solver

SolverType olb::FreeSurface::solver = SolverType::completePivotingLU

constexpr

Definition at line 63 of file freeSurfaceHelpers.h.