Skip to content

Commit

Permalink
Remove the global operations from SegregatedSolverBase. (idaholab#28819)
Browse files Browse the repository at this point in the history
  • Loading branch information
@grmnptr
grmnptr committed Oct 10, 2024
1 parent c473c8b commit 709e87d
Show file tree
Hide file tree
Showing 2 changed files with 1 addition and 270 deletions.
89 changes: 1 addition & 88 deletions modules/navier_stokes/include/executioners/SegregatedSolverBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
#include "INSFVRhieChowInterpolatorSegregated.h"
#include "PetscSupport.h"
#include "SolverParams.h"
#include "SegregatedSolverUtils.h"

#include "libmesh/petsc_vector.h"
#include "libmesh/petsc_matrix.h"
Expand Down Expand Up @@ -47,95 +48,7 @@ class SegregatedSolverBase : public Executioner
virtual void init() override;
virtual bool lastSolveConverged() const override { return _last_solve_converged; }

/**
* 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}
* }
* @param solution The solution vector
* @param mat The system matrix
* @param rhs The system right hand side
*/
Real computeNormalizationFactor(const NumericVector<Number> & solution,
const SparseMatrix<Number> & mat,
const NumericVector<Number> & rhs);

protected:
/**
* 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.
*
* @param matrix_in The matrix that needs to be relaxed
* @param relaxation_parameter The scale which described how much the matrix is relaxed
* @param diff_diagonal A vector holding the $A_{diag}-A_{diag, relaxed}$ entries for further
* use in the relaxation of the right hand side
*/
void relaxMatrix(SparseMatrix<Number> & matrix_in,
const Real relaxation_parameter,
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.
*
* @param rhs_in The unrelaxed right hand side that needs to be relaxed
* @param solution_in The solution
* @param diff_diagonal The diagonal correction used for the corresponding system matrix
*/
void relaxRightHandSide(NumericVector<Number> & rhs_in,
const NumericVector<Number> & solution_in,
const NumericVector<Number> & diff_diagonal);

/**
* Relax the update on a solution field using the following approach:
* $u = u_{old}+\lambda (u - u_{old})$
*
* @param vector_new The new solution vector
* @param vec_old The old solution vector
* @param relaxation_factor The lambda parameter in the expression above
*/
void relaxSolutionUpdate(NumericVector<Number> & vec_new,
const NumericVector<Number> & vec_old,
const Real relaxation_factor);

/**
* Limit a solution to its minimum and maximum bounds:
* $u = min(max(u, min_limit), max_limit)$
*
* @param system_in The system whose solution shall be limited
* @param min_limit = 0.0 The minimum limit for the solution
* @param max_limit = 1e10 The maximum limit for the solution
*/
void limitSolutionUpdate(NumericVector<Number> & solution,
const Real min_limit = std::numeric_limits<Real>::epsilon(),
const Real max_limit = 1e10);

/**
* 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.
* @param mx The matrix of the system which needs to be constrained
* @param rhs The right hand side of the system which needs to be constrained
* @param value The desired value for the solution field at a dof
* @param dof_id The ID of the dof which needs to be constrained
*/
void constrainSystem(SparseMatrix<Number> & mx,
NumericVector<Number> & rhs,
const Real desired_value,
const dof_id_type dof_id);

/**
* Find the ID of the degree of freedom which corresponds to the variable and
* a given point on the mesh
Expand Down
182 changes: 0 additions & 182 deletions modules/navier_stokes/src/executioners/SegregatedSolverBase.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -631,155 +631,6 @@ SegregatedSolverBase::init()
_pressure_pin_dof = findDoFID("pressure", getParam<Point>("pressure_pin_point"));
}

void
SegregatedSolverBase::relaxMatrix(SparseMatrix<Number> & matrix,
const Real relaxation_parameter,
NumericVector<Number> & diff_diagonal)
{
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
LIBMESH_CHKERR(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);
}

void
SegregatedSolverBase::relaxRightHandSide(NumericVector<Number> & rhs,
const NumericVector<Number> & solution,
const NumericVector<Number> & diff_diagonal)
{

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

Real
SegregatedSolverBase::computeNormalizationFactor(const NumericVector<Number> & solution,
const SparseMatrix<Number> & mat,
const NumericVector<Number> & rhs)
{
// 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);
}

void
SegregatedSolverBase::constrainSystem(SparseMatrix<Number> & mx,
NumericVector<Number> & rhs,
const Real desired_value,
const dof_id_type dof_id)
{
// 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);
}
}

dof_id_type
SegregatedSolverBase::findDoFID(const VariableName & var_name, const Point & point)
{
Expand Down Expand Up @@ -813,39 +664,6 @@ SegregatedSolverBase::findDoFID(const VariableName & var_name, const Point & poi
: DofObject::invalid_id;
}

void
SegregatedSolverBase::relaxSolutionUpdate(NumericVector<Number> & vec_new,
const NumericVector<Number> & vec_old,
const Real relaxation_factor)
{
// 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();

if (_print_fields)
{
_console << "Pressure solution" << std::endl;
vec_new.print();
_console << "Pressure solution old" << std::endl;
vec_old.print();
}
}

void
SegregatedSolverBase::limitSolutionUpdate(NumericVector<Number> & solution,
const Real min_limit,
const Real max_limit)
{
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();
}

bool
SegregatedSolverBase::converged(
const std::vector<std::pair<unsigned int, Real>> & its_and_residuals,
Expand Down

0 comments on commit 709e87d

Please sign in to comment.