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 213825 : relaxMatrix(SparseMatrix<Number> & matrix,
26 : const Real relaxation_parameter,
27 : NumericVector<Number> & diff_diagonal)
28 : {
29 213825 : PetscMatrix<Number> * mat = dynamic_cast<PetscMatrix<Number> *>(&matrix);
30 : mooseAssert(mat, "This should be a PetscMatrix!");
31 213825 : 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 213825 : *diff_diag = 0;
36 :
37 : // Get the diagonal of the matrix
38 213825 : mat->get_diagonal(*diff_diag);
39 :
40 : // Create a copy of the diagonal for later use and cast it
41 213825 : 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 213825 : 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 213825 : const unsigned int local_size = matrix.row_stop() - matrix.row_start();
60 213825 : std::vector<dof_id_type> indices(local_size);
61 213825 : std::iota(indices.begin(), indices.end(), matrix.row_start());
62 213825 : std::vector<Real> new_diagonal(local_size, 0.0);
63 :
64 : {
65 213825 : PetscVectorReader diff_diga_reader(*diff_diag);
66 17875728 : for (const auto row_i : make_range(local_size))
67 : {
68 17661903 : const unsigned int global_index = matrix.row_start() + row_i;
69 : std::vector<numeric_index_type> indices;
70 : std::vector<Real> values;
71 17661903 : mat->get_row(global_index, indices, values);
72 : const Real abs_sum = std::accumulate(
73 112602799 : values.cbegin(), values.cend(), 0.0, [](Real a, Real b) { return a + std::abs(b); });
74 17661903 : const Real abs_diagonal = std::abs(diff_diga_reader(global_index));
75 28243681 : new_diagonal[row_i] = inverse_relaxation * std::max(abs_sum - abs_diagonal, abs_diagonal);
76 : }
77 213825 : }
78 213825 : diff_diag->insert(new_diagonal, indices);
79 :
80 : // Time to modify the diagonal of the matrix. TODO: add this function to libmesh
81 213825 : LibmeshPetscCallA(mat->comm().get(), MatDiagonalSet(mat->mat(), diff_diag->vec(), INSERT_VALUES));
82 213825 : mat->close();
83 213825 : 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 213825 : diff_diag->add(-1.0, *original_diagonal);
88 213825 : }
89 :
90 : void
91 213825 : 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 213825 : auto working_vector = diff_diagonal.clone();
99 213825 : 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 213825 : rhs.add(*working_vector);
108 213825 : rhs.close();
109 213825 : }
110 :
111 : void
112 127385 : 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 127385 : vec_new.scale(relaxation_factor);
118 127385 : vec_new.add(1 - relaxation_factor, vec_old);
119 127385 : vec_new.close();
120 127385 : }
121 :
122 : void
123 4160 : limitSolutionUpdate(NumericVector<Number> & solution, const Real min_limit, const Real max_limit)
124 : {
125 4160 : PetscVector<Number> & solution_vector = dynamic_cast<PetscVector<Number> &>(solution);
126 : auto value = solution_vector.get_array();
127 :
128 352960 : for (auto i : make_range(solution_vector.local_size()))
129 697600 : value[i] = std::min(std::max(min_limit, value[i]), max_limit);
130 :
131 : solution_vector.restore_array();
132 4160 : }
133 :
134 : Real
135 337642 : 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 337642 : auto A_times_average_solution = solution.zero_clone();
153 337642 : 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 337642 : *A_times_solution = solution.sum() / solution.size();
159 337642 : mat.vector_mult(*A_times_average_solution, *A_times_solution);
160 337642 : mat.vector_mult(*A_times_solution, solution);
161 :
162 : // We create Ax-Ax_avg
163 337642 : 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 337642 : A_times_average_solution->add(-1.0, rhs);
166 337642 : A_times_solution->abs();
167 337642 : A_times_average_solution->abs();
168 :
169 : // Create |Ax-Ax_avg|+|b-Ax_avg|
170 337642 : 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 675284 : return (A_times_average_solution->l2_norm() + libMesh::TOLERANCE * libMesh::TOLERANCE);
177 337642 : }
178 :
179 : void
180 14311 : 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 14311 : if (dof_id >= mx.row_start() && dof_id < mx.row_stop())
189 : {
190 8210 : Real diag = mx(dof_id, dof_id);
191 8210 : rhs.add(dof_id, desired_value * diag);
192 8210 : mx.add(dof_id, dof_id, diag);
193 : }
194 14311 : }
195 :
196 : dof_id_type
197 426 : findPointDoFID(const MooseVariableFieldBase & variable, const MooseMesh & mesh, const Point & point)
198 : {
199 : // Find the element containing the point
200 426 : auto point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
201 426 : point_locator->enable_out_of_mesh_mode();
202 :
203 426 : 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 426 : variable.blockIDs().find(Moose::ANY_BLOCK_ID) == variable.blockIDs().end();
208 : const Elem * elem =
209 426 : 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 426 : const dof_id_type elem_id = elem ? elem->id() : DofObject::invalid_id;
214 426 : dof_id_type min_elem_id = elem_id;
215 426 : variable.sys().comm().min(min_elem_id);
216 :
217 426 : 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 426 : return min_elem_id == elem_id ? elem->dof_number(variable.sys().number(), var_num, 0)
225 426 : : DofObject::invalid_id;
226 426 : }
227 :
228 : bool
229 118344 : 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 133334 : for (const auto system_i : index_range(its_and_residuals))
239 : {
240 132521 : 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
|
