Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : // Moose includes
11 : #include "SegregatedSolverUtils.h"
12 : #include "PetscVectorReader.h"
13 : #include "MooseVariableFieldBase.h"
14 : #include "SystemBase.h"
15 :
16 : // Libmesh includes
17 : #include "libmesh/enum_point_locator_type.h"
18 : #include "libmesh/system.h"
19 :
20 : namespace NS
21 : {
22 : namespace FV
23 : {
24 : void
25 494404 : relaxMatrix(SparseMatrix<Number> & matrix,
26 : const Real relaxation_parameter,
27 : NumericVector<Number> & diff_diagonal)
28 : {
29 494404 : PetscMatrix<Number> * mat = dynamic_cast<PetscMatrix<Number> *>(&matrix);
30 : mooseAssert(mat, "This should be a PetscMatrix!");
31 494404 : PetscVector<Number> * diff_diag = dynamic_cast<PetscVector<Number> *>(&diff_diagonal);
32 : mooseAssert(diff_diag, "This should be a PetscVector!");
33 :
34 : // Zero the diagonal difference vector
35 494404 : *diff_diag = 0;
36 :
37 : // Get the diagonal of the matrix
38 494404 : mat->get_diagonal(*diff_diag);
39 :
40 : // Create a copy of the diagonal for later use and cast it
41 494404 : auto original_diagonal = diff_diag->clone();
42 :
43 : // We cache the inverse of the relaxation parameter because doing divisions might
44 : // be more expensive for every row
45 494404 : const Real inverse_relaxation = 1 / relaxation_parameter;
46 :
47 : // Now we loop over the matrix row by row and sum the absolute values of the
48 : // offdiagonal values, If these sums are larger than the diagonal entry,
49 : // we switch the diagonal value with the sum. At the same time we increase the
50 : // diagonal by dividing it with the relaxation parameter. So the new diagonal will be:
51 : // D* = 1/lambda*max(|D|,sum(|Offdiag|))
52 : // For more information see
53 : //
54 : // Juretic, Franjo. Error analysis in finite volume CFD. Diss.
55 : // Imperial College London (University of London), 2005.
56 : //
57 : // The trickery comes with storing everything in the diff-diagonal vector
58 : // to avoid the allocation and manipulation of a third vector
59 494404 : const unsigned int local_size = matrix.row_stop() - matrix.row_start();
60 494404 : std::vector<dof_id_type> indices(local_size);
61 494404 : std::iota(indices.begin(), indices.end(), matrix.row_start());
62 494404 : std::vector<Real> new_diagonal(local_size, 0.0);
63 :
64 : {
65 494404 : PetscVectorReader diff_diga_reader(*diff_diag);
66 37232190 : for (const auto row_i : make_range(local_size))
67 : {
68 36737786 : const unsigned int global_index = matrix.row_start() + row_i;
69 : std::vector<numeric_index_type> indices;
70 : std::vector<Real> values;
71 36737786 : mat->get_row(global_index, indices, values);
72 : const Real abs_sum = std::accumulate(
73 212313996 : values.cbegin(), values.cend(), 0.0, [](Real a, Real b) { return a + std::abs(b); });
74 36737786 : const Real abs_diagonal = std::abs(diff_diga_reader(global_index));
75 58003326 : new_diagonal[row_i] = inverse_relaxation * std::max(abs_sum - abs_diagonal, abs_diagonal);
76 36737786 : }
77 494404 : }
78 494404 : diff_diag->insert(new_diagonal, indices);
79 :
80 : // Time to modify the diagonal of the matrix. TODO: add this function to libmesh
81 494404 : LibmeshPetscCallA(mat->comm().get(), MatDiagonalSet(mat->mat(), diff_diag->vec(), INSERT_VALUES));
82 494404 : mat->close();
83 494404 : diff_diag->close();
84 :
85 : // Finally, we can create (D*-D) vector which is used for the relaxation of the
86 : // right hand side later
87 494404 : diff_diag->add(-1.0, *original_diagonal);
88 494404 : }
89 :
90 : void
91 494404 : relaxRightHandSide(NumericVector<Number> & rhs,
92 : const NumericVector<Number> & solution,
93 : const NumericVector<Number> & diff_diagonal)
94 : {
95 :
96 : // We need a working vector here to make sure we don't modify the
97 : // (D*-D) vector
98 494404 : auto working_vector = diff_diagonal.clone();
99 494404 : working_vector->pointwise_mult(solution, *working_vector);
100 :
101 : // The correction to the right hand side is just
102 : // (D*-D)*old_solution
103 : // For more information see
104 : //
105 : // Juretic, Franjo. Error analysis in finite volume CFD. Diss.
106 : // Imperial College London (University of London), 2005.
107 494404 : rhs.add(*working_vector);
108 494404 : rhs.close();
109 494404 : }
110 :
111 : void
112 257792 : relaxSolutionUpdate(NumericVector<Number> & vec_new,
113 : const NumericVector<Number> & vec_old,
114 : const Real relaxation_factor)
115 : {
116 : // The relaxation is just u = lambda * u* + (1-lambda) u_old
117 257792 : vec_new.scale(relaxation_factor);
118 257792 : vec_new.add(1 - relaxation_factor, vec_old);
119 257792 : vec_new.close();
120 257792 : }
121 :
122 : void
123 55628 : limitSolutionUpdate(NumericVector<Number> & solution, const Real min_limit, const Real max_limit)
124 : {
125 55628 : PetscVector<Number> & solution_vector = dynamic_cast<PetscVector<Number> &>(solution);
126 : auto value = solution_vector.get_array();
127 :
128 1734172 : for (auto i : make_range(solution_vector.local_size()))
129 3357088 : value[i] = std::min(std::max(min_limit, value[i]), max_limit);
130 :
131 : solution_vector.restore_array();
132 55628 : }
133 :
134 : Real
135 727508 : computeNormalizationFactor(const NumericVector<Number> & solution,
136 : const SparseMatrix<Number> & mat,
137 : const NumericVector<Number> & rhs)
138 : {
139 : // This function is based on the description provided here:
140 : // @article{greenshields2022notes,
141 : // title={Notes on computational fluid dynamics: General principles},
142 : // author={Greenshields, Christopher J and Weller, Henry G},
143 : // journal={(No Title)},
144 : // year={2022}
145 : // }
146 : // so basically we normalize the residual with the following number:
147 : // sum(|Ax-Ax_avg|+|b-Ax_avg|)
148 : // where A is the system matrix, b is the system right hand side while x and x_avg are
149 : // the solution and average solution vectors
150 :
151 : // We create a vector for Ax_avg and Ax
152 727508 : auto A_times_average_solution = solution.zero_clone();
153 727508 : auto A_times_solution = solution.zero_clone();
154 :
155 : // Beware, trickery here! To avoid allocating unused vectors, we
156 : // first compute Ax_avg using the storage used for Ax, then we
157 : // overwrite Ax with the right value
158 727508 : *A_times_solution = solution.sum() / solution.size();
159 727508 : mat.vector_mult(*A_times_average_solution, *A_times_solution);
160 727508 : mat.vector_mult(*A_times_solution, solution);
161 :
162 : // We create Ax-Ax_avg
163 727508 : A_times_solution->add(-1.0, *A_times_average_solution);
164 : // We create Ax_avg - b (ordering shouldn't matter we will take absolute value soon)
165 727508 : A_times_average_solution->add(-1.0, rhs);
166 727508 : A_times_solution->abs();
167 727508 : A_times_average_solution->abs();
168 :
169 : // Create |Ax-Ax_avg|+|b-Ax_avg|
170 727508 : A_times_average_solution->add(*A_times_solution);
171 :
172 : // Since use the l2 norm of the solution vectors in the linear solver, we will
173 : // make this consistent and use the l2 norm of the vector. We add a small number to
174 : // avoid normalizing with 0.
175 : // TODO: Would be nice to see if we can do l1 norms in the linear solve.
176 1455016 : return (A_times_average_solution->l2_norm() + libMesh::TOLERANCE * libMesh::TOLERANCE);
177 727508 : }
178 :
179 : void
180 67765 : constrainSystem(SparseMatrix<Number> & mx,
181 : NumericVector<Number> & rhs,
182 : const Real desired_value,
183 : const dof_id_type dof_id)
184 : {
185 : // Modify the given matrix and right hand side. We use the matrix diagonal
186 : // to enforce the constraint instead of 1, to make sure we don't mess up the matrix conditioning
187 : // too much.
188 67765 : if (dof_id >= mx.row_start() && dof_id < mx.row_stop())
189 : {
190 30617 : Real diag = mx(dof_id, dof_id);
191 30617 : rhs.add(dof_id, desired_value * diag);
192 30617 : mx.add(dof_id, dof_id, diag);
193 : }
194 67765 : }
195 :
196 : dof_id_type
197 830 : findPointDoFID(const MooseVariableFieldBase & variable, const MooseMesh & mesh, const Point & point)
198 : {
199 : // Find the element containing the point
200 830 : auto point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
201 830 : point_locator->enable_out_of_mesh_mode();
202 :
203 830 : unsigned int var_num = variable.sys().system().variable_number(variable.name());
204 :
205 : // We only check in the restricted blocks, if needed
206 : const bool block_restricted =
207 830 : variable.blockIDs().find(Moose::ANY_BLOCK_ID) == variable.blockIDs().end();
208 : const Elem * elem =
209 830 : block_restricted ? (*point_locator)(point, &variable.blockIDs()) : (*point_locator)(point);
210 :
211 : // We communicate the results and if there is conflict between processes,
212 : // the minimum cell ID is chosen
213 830 : const dof_id_type elem_id = elem ? elem->id() : DofObject::invalid_id;
214 830 : dof_id_type min_elem_id = elem_id;
215 830 : variable.sys().comm().min(min_elem_id);
216 :
217 830 : if (min_elem_id == DofObject::invalid_id)
218 0 : mooseError("Variable ",
219 0 : variable.name(),
220 : " is not defined at ",
221 0 : Moose::stringify(point),
222 : "! Try alleviating block restrictions or using another point!");
223 :
224 830 : return min_elem_id == elem_id ? elem->dof_number(variable.sys().number(), var_num, 0)
225 830 : : DofObject::invalid_id;
226 830 : }
227 :
228 : bool
229 218647 : converged(const std::vector<std::pair<unsigned int, Real>> & its_and_residuals,
230 : const std::vector<Real> & abs_tolerances)
231 : {
232 : mooseAssert(its_and_residuals.size() == abs_tolerances.size(),
233 : "The number of residuals should (now " + std::to_string(its_and_residuals.size()) +
234 : ") be the same as the number of tolerances (" +
235 : std::to_string(abs_tolerances.size()) + ")!");
236 :
237 : bool converged = true;
238 257653 : for (const auto system_i : index_range(its_and_residuals))
239 : {
240 256093 : converged = converged && (its_and_residuals[system_i].second < abs_tolerances[system_i]);
241 : if (!converged)
242 : return converged;
243 : }
244 : return converged;
245 : }
246 :
247 : } // End FV namespace
248 : } // End Moose namespace
|
