NS::FV Namespace Reference

Functions
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. More...
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. More...
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}+ (u - u_{old})$. More...
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)$. More...
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. More...
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. More...
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. More...
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. More...

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: {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.

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

{
  // 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);
}

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.

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

{
  // 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);
  }
}

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.

Referenced by MechanicalContactConstraint::AugmentedLagrangianContactConverged(), GeneralizedReturnMappingSolutionTempl< is_ad >::convergedAcceptable(), SingleVariableReturnMappingSolutionTempl< is_ad >::convergedAcceptable(), SingularShapeTensorEliminatorUserObjectPD::execute(), externalPETScDiffusionFDMSolve(), GeneralizedReturnMappingSolutionTempl< is_ad >::internalSolve(), SingleVariableReturnMappingSolutionTempl< is_ad >::internalSolve(), FluidPropertiesUtils::NewtonSolve(), FluidPropertiesUtils::NewtonSolve2D(), ExplicitMixedOrder::performExplicitSolve(), LinearAssemblySegregatedSolve::solve(), SIMPLESolveNonlinearAssembly::solve(), ExplicitMixedOrder::solve(), and AdjointTransientSolve::solve().

{
  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;
}

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.

Referenced by SIMPLESolveBase::setupPressurePin().

{
  // 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;
}

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.

Referenced by SIMPLESolveNonlinearAssembly::solve().

{
  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();
}

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.

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

{
  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);
}

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.

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

{

  // 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();
}

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}+ (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.

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

{
  // 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();
}