NS::FV Namespace Reference

Classes
class	CHTHandler
	This class provides an interface for managing conjugate heat transfer (CHT) between fluid and solid domains. More...
struct	PressureJumpInterfaceData
	Complete face operator obtained by eliminating the two one-sided interface pressures. More...

Functions
PressureJumpInterfaceData	pressureJumpInterfaceData (const RealVectorValue &normal, const Point &elem_to_face, const Point &face_to_neighbor, const RealVectorValue &elem_diffusion, const RealVectorValue &neighbor_diffusion, const RealVectorValue &elem_gradient, const RealVectorValue &neighbor_gradient, const RealVectorValue &elem_rho_hbya, const RealVectorValue &neighbor_rho_hbya, const Real face_area, const bool use_nonorthogonal_correction, const Real elem_jump)
	Construct the conservative half-cell interface coefficients for a diagonal diffusion tensor.
Real	pressureJumpHbyAFlux (const PressureJumpInterfaceData &data)
	Compute the H/A part of the face mass flux after eliminating the interface pressures.
Real	pressureJumpFlux (const PressureJumpInterfaceData &data, const Real elem_pressure, const Real neighbor_pressure)
	Compute the pressure part of the element-to-neighbor flux after eliminating the interface pressures.
Real	pressureJumpOneSidedFaceValue (const PressureJumpInterfaceData &data, const Real cell_pressure, const Real mass_flux, const bool elem_side)
	Recover one side's face pressure from the common oriented face mass flux.
void	relaxMatrix (SparseMatrix< Number > &matrix_in, const Real relaxation_parameter, NumericVector< Number > &diff_diagonal)
	Relax the matrix to ensure diagonal dominance, we hold onto the difference in diagonals for later use in relaxing the right hand side.
void	relaxRightHandSide (NumericVector< Number > &rhs_in, const NumericVector< Number > &solution_in, const NumericVector< Number > &diff_diagonal)
	Relax the right hand side of an equation, this needs to be called once and the system matrix has been relaxed and the field describing the difference in the diagonals of the system matrix is already available.
void	relaxSolutionUpdate (NumericVector< Number > &vec_new, const NumericVector< Number > &vec_old, const Real relaxation_factor)
	Relax the update on a solution field using the following approach: $u = u_{old}+\lambda (u - u_{old})$.
void	limitSolutionUpdate (NumericVector< Number > &solution, const Real min_limit=std::numeric_limits< Real >::epsilon(), const Real max_limit=1e10)
	Limit a solution to its minimum and maximum bounds: $u = min(max(u, min_limit), max_limit)$.
Real	computeNormalizationFactor (const NumericVector< Number > &solution, const SparseMatrix< Number > &mat, const NumericVector< Number > &rhs)
	Compute a normalization factor which is applied to the linear residual to determine convergence.
void	constrainSystem (SparseMatrix< Number > &mx, NumericVector< Number > &rhs, const Real desired_value, const dof_id_type dof_id)
	Implicitly constrain the system by adding a factor*(u-u_desired) to it at a desired dof value.
dof_id_type	findPointDoFID (const MooseVariableFieldBase &variable, const MooseMesh &mesh, const Point &point)
	Find the ID of the degree of freedom which corresponds to the variable and a given point on the mesh.
bool	converged (const std::vector< std::pair< unsigned int, Real > > &residuals, const std::vector< Real > &abs_tolerances)
	Based on the residuals, determine if the iterative process converged or not.

Function Documentation

computeNormalizationFactor()

Real NS::FV::computeNormalizationFactor	(	const NumericVector< Number > &	solution,
		const SparseMatrix< Number > &	mat,
		const NumericVector< Number > &	rhs
	)

Compute a normalization factor which is applied to the linear residual to determine convergence.

This function is based on the description provided here: @article{greenshields2022notes, title={Notes on computational fluid dynamics: General principles}, author={Greenshields, Christopher J and Weller, Henry G}, journal={(No Title)}, year={2022} }

Parameters

solution	The solution vector
mat	The system matrix
rhs	The system right hand side

Definition at line 135 of file SegregatedSolverUtils.C.

{
  // This function is based on the description provided here:
  // @article{greenshields2022notes,
  // title={Notes on computational fluid dynamics: General principles},
  // author={Greenshields, Christopher J and Weller, Henry G},
  // journal={(No Title)},
  // year={2022}
  // }
  // so basically we normalize the residual with the following number:
  // sum(|Ax-Ax_avg|+|b-Ax_avg|)
  // where A is the system matrix, b is the system right hand side while x and x_avg are
  // the solution and average solution vectors
 
  // We create a vector for Ax_avg and Ax
  auto A_times_average_solution = solution.zero_clone();
  auto A_times_solution = solution.zero_clone();
 
  // Beware, trickery here! To avoid allocating unused vectors, we
  // first compute Ax_avg using the storage used for Ax, then we
  // overwrite Ax with the right value
  *A_times_solution = solution.sum() / solution.size();
  mat.vector_mult(*A_times_average_solution, *A_times_solution);
  mat.vector_mult(*A_times_solution, solution);
 
  // We create Ax-Ax_avg
  A_times_solution->add(-1.0, *A_times_average_solution);
  // We create Ax_avg - b (ordering shouldn't matter we will take absolute value soon)
  A_times_average_solution->add(-1.0, rhs);
  A_times_solution->abs();
  A_times_average_solution->abs();
 
  // Create |Ax-Ax_avg|+|b-Ax_avg|
  A_times_average_solution->add(*A_times_solution);
 
  // Since use the l2 norm of the solution vectors in the linear solver, we will
  // make this consistent and use the l2 norm of the vector. We add a small number to
  // avoid normalizing with 0.
  // TODO: Would be nice to see if we can do l1 norms in the linear solve.
  return (A_times_average_solution->l2_norm() + libMesh::TOLERANCE * libMesh::TOLERANCE);
}

Referenced by LinearAssemblySegregatedSolve::solveAdvectedSystem(), SIMPLESolveNonlinearAssembly::solveAdvectedSystem(), LinearAssemblySegregatedSolve::solveMomentumPredictor(), SIMPLESolveNonlinearAssembly::solveMomentumPredictor(), LinearAssemblySegregatedSolve::solvePressureCorrector(), SIMPLESolveNonlinearAssembly::solvePressureCorrector(), LinearAssemblySegregatedSolve::solveSolidEnergy(), and SIMPLESolveNonlinearAssembly::solveSolidEnergySystem().

constrainSystem()

void NS::FV::constrainSystem	(	SparseMatrix< Number > &	mx,
		NumericVector< Number > &	rhs,
		const Real	desired_value,
		const dof_id_type	dof_id
	)

Implicitly constrain the system by adding a factor*(u-u_desired) to it at a desired dof value.

To make sure the conditioning of the matrix does not change significantly, factor is chosen to be the diagonal component of the matrix coefficients for a given dof.

Parameters

mx	The matrix of the system which needs to be constrained
rhs	The right hand side of the system which needs to be constrained
value	The desired value for the solution field at a dof
dof_id	The ID of the dof which needs to be constrained

Definition at line 180 of file SegregatedSolverUtils.C.

{
  // Modify the given matrix and right hand side. We use the matrix diagonal
  // to enforce the constraint instead of 1, to make sure we don't mess up the matrix conditioning
  // too much.
  if (dof_id >= mx.row_start() && dof_id < mx.row_stop())
  {
    Real diag = mx(dof_id, dof_id);
    rhs.add(dof_id, desired_value * diag);
    mx.add(dof_id, dof_id, diag);
  }
}

Referenced by LinearAssemblySegregatedSolve::solvePressureCorrector(), and SIMPLESolveNonlinearAssembly::solvePressureCorrector().

converged()

bool NS::FV::converged	(	const std::vector< std::pair< unsigned int, Real > > &	residuals,
		const std::vector< Real > &	abs_tolerances
	)

Based on the residuals, determine if the iterative process converged or not.

Parameters

residuals	The current (number of iterations, residual) pairs
abs_tolerances	The corresponding absolute tolerances.

Definition at line 229 of file SegregatedSolverUtils.C.

{
  mooseAssert(its_and_residuals.size() == abs_tolerances.size(),
              "The number of residuals should (now " + std::to_string(its_and_residuals.size()) +
                  ") be the same as the number of tolerances (" +
                  std::to_string(abs_tolerances.size()) + ")!");
 
  bool converged = true;
  for (const auto system_i : index_range(its_and_residuals))
  {
    converged = converged && (its_and_residuals[system_i].second < abs_tolerances[system_i]);
    if (!converged)
      return converged;
  }
  return converged;
}

Referenced by converged(), LinearAssemblySegregatedSolve::solve(), and SIMPLESolveNonlinearAssembly::solve().

findPointDoFID()

dof_id_type NS::FV::findPointDoFID	(	const MooseVariableFieldBase &	variable,
		const MooseMesh &	mesh,
		const Point &	point
	)

Find the ID of the degree of freedom which corresponds to the variable and a given point on the mesh.

Parameters

variable	Reference to the moose variable whose dof should be fetched
mesh	The moose mesh where the element is
point	The point on the mesh

Definition at line 197 of file SegregatedSolverUtils.C.

{
  // Find the element containing the point
  auto point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
  point_locator->enable_out_of_mesh_mode();
 
  unsigned int var_num = variable.sys().system().variable_number(variable.name());
 
  // We only check in the restricted blocks, if needed
  const bool block_restricted =
      variable.blockIDs().find(Moose::ANY_BLOCK_ID) == variable.blockIDs().end();
  const Elem * elem =
      block_restricted ? (*point_locator)(point, &variable.blockIDs()) : (*point_locator)(point);
 
  // We communicate the results and if there is conflict between processes,
  // the minimum cell ID is chosen
  const dof_id_type elem_id = elem ? elem->id() : DofObject::invalid_id;
  dof_id_type min_elem_id = elem_id;
  variable.sys().comm().min(min_elem_id);
 
  if (min_elem_id == DofObject::invalid_id)
    mooseError("Variable ",
               variable.name(),
               " is not defined at ",
               Moose::stringify(point),
               "! Try alleviating block restrictions or using another point!");
 
  return min_elem_id == elem_id ? elem->dof_number(variable.sys().number(), var_num, 0)
                                : DofObject::invalid_id;
}

Referenced by SIMPLESolveBase::setupPressurePin().

limitSolutionUpdate()

void NS::FV::limitSolutionUpdate	(	NumericVector< Number > &	solution,
		const Real	min_limit = std::numeric_limits<Real>::epsilon(),
		const Real	max_limit = 1e10
	)

Limit a solution to its minimum and maximum bounds: $u = min(max(u, min_limit), max_limit)$.

Parameters

system_in	The system whose solution shall be limited
min_limit	= 0.0 The minimum limit for the solution
max_limit	= 1e10 The maximum limit for the solution

Definition at line 123 of file SegregatedSolverUtils.C.

{
  PetscVector<Number> & solution_vector = dynamic_cast<PetscVector<Number> &>(solution);
  auto value = solution_vector.get_array();
 
  for (auto i : make_range(solution_vector.local_size()))
    value[i] = std::min(std::max(min_limit, value[i]), max_limit);
 
  solution_vector.restore_array();
}

Referenced by SIMPLESolveNonlinearAssembly::solve(), and LinearAssemblySegregatedSolve::solveAdvectedSystem().

pressureJumpFlux()

Real NS::FV::pressureJumpFlux ( const PressureJumpInterfaceData &  data, const Real  elem_pressure, const Real  neighbor_pressure  )

inline

Compute the pressure part of the element-to-neighbor flux after eliminating the interface pressures.

The face mass flux is this minus pressureJumpHbyAFlux.

Definition at line 138 of file PressureJumpInterface.h.

{
  return data.transmissibility * (elem_pressure - neighbor_pressure - data.elem_jump) -
         data.correction;
}

Referenced by FVPressureJumpGreenGaussGradient::internalFaceValues(), TEST(), TEST(), TEST(), and TEST().

pressureJumpHbyAFlux()

Real NS::FV::pressureJumpHbyAFlux ( const PressureJumpInterfaceData &  data )

inline

Compute the H/A part of the face mass flux after eliminating the interface pressures.

Each side carries mdot = -h_s + (pressure flux of side s), so continuity of mdot weights the side H/A fluxes h_s with the same series weights T / T_s as the pressure.

Definition at line 129 of file PressureJumpInterface.h.

{
  return data.transmissibility * (data.elem_hbya_flux / data.elem_conductance +
                                  data.neighbor_hbya_flux / data.neighbor_conductance);
}

Referenced by PorousRhieChowMassFlux::computeHbyA(), FVPressureJumpGreenGaussGradient::internalFaceValues(), and TEST().

pressureJumpInterfaceData()

PressureJumpInterfaceData NS::FV::pressureJumpInterfaceData ( const RealVectorValue &  normal, const Point &  elem_to_face, const Point &  face_to_neighbor, const RealVectorValue &  elem_diffusion, const RealVectorValue &  neighbor_diffusion, const RealVectorValue &  elem_gradient, const RealVectorValue &  neighbor_gradient, const RealVectorValue &  elem_rho_hbya, const RealVectorValue &  neighbor_rho_hbya, const Real  face_area, const bool  use_nonorthogonal_correction, const Real  elem_jump  )

inline

Construct the conservative half-cell interface coefficients for a diagonal diffusion tensor.

Both half-cell direction vectors point in the element-to-neighbor orientation.

Definition at line 48 of file PressureJumpInterface.h.

{
  PressureJumpInterfaceData data;
  data.is_baffle = true;
  data.elem_jump = elem_jump;
  data.face_area = face_area;
 
  const auto side_data = [&](const Point & half_cell,
                             const RealVectorValue & diffusion,
                             const RealVectorValue & gradient,
                             Real & conductance,
                             Real & correction,
                             Real & normal_distance,
                             Real & normal_diffusion)
  {
    RealVectorValue diffusion_normal;
    for (const auto i : make_range(Moose::dim))
      diffusion_normal(i) = diffusion(i) * normal(i);
 
    normal_diffusion = normal * diffusion_normal;
    normal_distance = half_cell * normal;
    if (!std::isfinite(normal_distance) || normal_distance <= 0.0 ||
        !std::isfinite(normal_diffusion) || normal_diffusion <= 0.0 || face_area <= 0.0)
      return false;
 
    conductance = face_area * normal_diffusion / normal_distance;
 
    RealVectorValue correction_vector = diffusion_normal - normal_diffusion * normal;
    if (use_nonorthogonal_correction)
      correction_vector += normal_diffusion * (normal - half_cell / normal_distance);
 
    correction = face_area * correction_vector * gradient;
    return std::isfinite(conductance) && conductance > 0.0 && std::isfinite(correction);
  };
 
  const bool elem_valid = side_data(elem_to_face,
                                    elem_diffusion,
                                    elem_gradient,
                                    data.elem_conductance,
                                    data.elem_correction,
                                    data.elem_normal_distance,
                                    data.elem_normal_diffusion);
  const bool neighbor_valid = side_data(face_to_neighbor,
                                        neighbor_diffusion,
                                        neighbor_gradient,
                                        data.neighbor_conductance,
                                        data.neighbor_correction,
                                        data.neighbor_normal_distance,
                                        data.neighbor_normal_diffusion);
  if (!elem_valid || !neighbor_valid)
    return data;
 
  data.elem_hbya_flux = face_area * (elem_rho_hbya * normal);
  data.neighbor_hbya_flux = face_area * (neighbor_rho_hbya * normal);
  // conductance in series
  data.transmissibility = 1.0 / (1.0 / data.elem_conductance + 1.0 / data.neighbor_conductance);
  data.correction = data.transmissibility * (data.elem_correction / data.elem_conductance +
                                             data.neighbor_correction / data.neighbor_conductance);
  data.valid = std::isfinite(data.transmissibility) && data.transmissibility > 0.0 &&
               std::isfinite(data.correction) && std::isfinite(data.elem_hbya_flux) &&
               std::isfinite(data.neighbor_hbya_flux);
  return data;
}

Referenced by PorousRhieChowMassFlux::pressureJumpInterfaceData(), TEST(), TEST(), TEST(), TEST(), TEST(), and TEST().

pressureJumpOneSidedFaceValue()

Real NS::FV::pressureJumpOneSidedFaceValue ( const PressureJumpInterfaceData &  data, const Real  cell_pressure, const Real  mass_flux, const bool  elem_side  )

inline

Recover one side's face pressure from the common oriented face mass flux.

The face value of the pressure variable interpolates the two cell pressures, which is valid only where the pressure is continuous. On a baffle face, any consumer that needs the face pressure of one side, such as the surface sum of a Green-Gauss gradient, must use this value instead.

Definition at line 154 of file PressureJumpInterface.h.

{
  if (elem_side)
    return cell_pressure -
           (mass_flux + data.elem_hbya_flux + data.elem_correction) / data.elem_conductance;
 
  return cell_pressure + (mass_flux + data.neighbor_hbya_flux + data.neighbor_correction) /
                             data.neighbor_conductance;
}

Referenced by FVPressureJumpGreenGaussGradient::internalFaceValues(), TEST(), TEST(), and TEST().

relaxMatrix()

void NS::FV::relaxMatrix	(	SparseMatrix< Number > &	matrix_in,
		const Real	relaxation_parameter,
		NumericVector< Number > &	diff_diagonal
	)

Relax the matrix to ensure diagonal dominance, we hold onto the difference in diagonals for later use in relaxing the right hand side.

For the details of this relaxation process, see

Juretic, Franjo. Error analysis in finite volume CFD. Diss. Imperial College London (University of London), 2005.

Parameters

matrix_in	The matrix that needs to be relaxed
relaxation_parameter	The scale which described how much the matrix is relaxed
diff_diagonal	A vector holding the $A_{diag}-A_{diag, relaxed}$ entries for further use in the relaxation of the right hand side

Definition at line 25 of file SegregatedSolverUtils.C.

{
  PetscMatrix<Number> * mat = dynamic_cast<PetscMatrix<Number> *>(&matrix);
  mooseAssert(mat, "This should be a PetscMatrix!");
  PetscVector<Number> * diff_diag = dynamic_cast<PetscVector<Number> *>(&diff_diagonal);
  mooseAssert(diff_diag, "This should be a PetscVector!");
 
  // Zero the diagonal difference vector
  *diff_diag = 0;
 
  // Get the diagonal of the matrix
  mat->get_diagonal(*diff_diag);
 
  // Create a copy of the diagonal for later use and cast it
  auto original_diagonal = diff_diag->clone();
 
  // We cache the inverse of the relaxation parameter because doing divisions might
  // be more expensive for every row
  const Real inverse_relaxation = 1 / relaxation_parameter;
 
  // Now we loop over the matrix row by row and sum the absolute values of the
  // offdiagonal values, If these sums are larger than the diagonal entry,
  // we switch the diagonal value with the sum. At the same time we increase the
  // diagonal by dividing it with the relaxation parameter. So the new diagonal will be:
  // D* = 1/lambda*max(|D|,sum(|Offdiag|))
  // For more information see
  //
  // Juretic, Franjo. Error analysis in finite volume CFD. Diss.
  // Imperial College London (University of London), 2005.
  //
  // The trickery comes with storing everything in the diff-diagonal vector
  // to avoid the allocation and manipulation of a third vector
  const unsigned int local_size = matrix.row_stop() - matrix.row_start();
  std::vector<dof_id_type> indices(local_size);
  std::iota(indices.begin(), indices.end(), matrix.row_start());
  std::vector<Real> new_diagonal(local_size, 0.0);
 
  {
    PetscVectorReader diff_diga_reader(*diff_diag);
    for (const auto row_i : make_range(local_size))
    {
      const unsigned int global_index = matrix.row_start() + row_i;
      std::vector<numeric_index_type> indices;
      std::vector<Real> values;
      mat->get_row(global_index, indices, values);
      const Real abs_sum = std::accumulate(
          values.cbegin(), values.cend(), 0.0, [](Real a, Real b) { return a + std::abs(b); });
      const Real abs_diagonal = std::abs(diff_diga_reader(global_index));
      new_diagonal[row_i] = inverse_relaxation * std::max(abs_sum - abs_diagonal, abs_diagonal);
    }
  }
  diff_diag->insert(new_diagonal, indices);
 
  // Time to modify the diagonal of the matrix. TODO: add this function to libmesh
  LibmeshPetscCallA(mat->comm().get(), MatDiagonalSet(mat->mat(), diff_diag->vec(), INSERT_VALUES));
  mat->close();
  diff_diag->close();
 
  // Finally, we can create (D*-D) vector which is used for the relaxation of the
  // right hand side later
  diff_diag->add(-1.0, *original_diagonal);
}

Referenced by LinearAssemblySegregatedSolve::solveAdvectedSystem(), SIMPLESolveNonlinearAssembly::solveAdvectedSystem(), LinearAssemblySegregatedSolve::solveMomentumPredictor(), and SIMPLESolveNonlinearAssembly::solveMomentumPredictor().

relaxRightHandSide()

void NS::FV::relaxRightHandSide	(	NumericVector< Number > &	rhs_in,
		const NumericVector< Number > &	solution_in,
		const NumericVector< Number > &	diff_diagonal
	)

Relax the right hand side of an equation, this needs to be called once and the system matrix has been relaxed and the field describing the difference in the diagonals of the system matrix is already available.

The relaxation process needs modification to both the system matrix and the right hand side. For more information see:

Juretic, Franjo. Error analysis in finite volume CFD. Diss. Imperial College London (University of London), 2005.

Parameters

rhs_in	The unrelaxed right hand side that needs to be relaxed
solution_in	The solution
diff_diagonal	The diagonal correction used for the corresponding system matrix

Definition at line 91 of file SegregatedSolverUtils.C.

{
 
  // We need a working vector here to make sure we don't modify the
  // (D*-D) vector
  auto working_vector = diff_diagonal.clone();
  working_vector->pointwise_mult(solution, *working_vector);
 
  // The correction to the right hand side is just
  // (D*-D)*old_solution
  // For more information see
  //
  // Juretic, Franjo. Error analysis in finite volume CFD. Diss.
  // Imperial College London (University of London), 2005.
  rhs.add(*working_vector);
  rhs.close();
}

Referenced by LinearAssemblySegregatedSolve::solveAdvectedSystem(), SIMPLESolveNonlinearAssembly::solveAdvectedSystem(), LinearAssemblySegregatedSolve::solveMomentumPredictor(), and SIMPLESolveNonlinearAssembly::solveMomentumPredictor().

relaxSolutionUpdate()

void NS::FV::relaxSolutionUpdate	(	NumericVector< Number > &	vec_new,
		const NumericVector< Number > &	vec_old,
		const Real	relaxation_factor
	)

Relax the update on a solution field using the following approach: $u = u_{old}+\lambda (u - u_{old})$.

Parameters

vector_new	The new solution vector
vec_old	The old solution vector
relaxation_factor	The lambda parameter in the expression above

Definition at line 112 of file SegregatedSolverUtils.C.

{
  // The relaxation is just u = lambda * u* + (1-lambda) u_old
  vec_new.scale(relaxation_factor);
  vec_new.add(1 - relaxation_factor, vec_old);
  vec_new.close();
}

Referenced by LinearAssemblySegregatedSolve::correctVelocity(), SIMPLESolveNonlinearAssembly::solve(), and LinearAssemblySegregatedSolve::solveAdvectedSystem().