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 : #include "NavierStokesProblem.h"
11 : #include "NonlinearSystemBase.h"
12 : #include "libmesh/petsc_matrix.h"
13 :
14 : using namespace libMesh;
15 :
16 : registerMooseObject("NavierStokesApp", NavierStokesProblem);
17 :
18 : InputParameters
19 330 : NavierStokesProblem::validParams()
20 : {
21 330 : InputParameters params = FEProblem::validParams();
22 330 : params.addClassDescription("A problem that handles Schur complement preconditioning of the "
23 : "incompressible Navier-Stokes equations");
24 660 : params.addParam<TagName>(
25 : "mass_matrix", "", "The matrix tag name corresponding to the mass matrix.");
26 660 : params.addParam<TagName>(
27 : "L_matrix",
28 : "",
29 : "The matrix tag name corresponding to the diffusive part of the velocity equations.");
30 660 : params.addParam<std::vector<unsigned int>>(
31 : "schur_fs_index",
32 : {},
33 : "if not provided then the top field split is assumed to be the "
34 : "Schur split. This is a vector to allow recursive nesting");
35 660 : params.addParam<bool>("use_pressure_mass_matrix",
36 660 : false,
37 : "Whether to just use the pressure mass matrix as the preconditioner for "
38 : "the Schur complement");
39 660 : params.addParam<bool>(
40 : "commute_lsc",
41 660 : false,
42 : "Whether to use the commuted form of the LSC preconditioner, created by Olshanskii");
43 330 : return params;
44 0 : }
45 :
46 165 : NavierStokesProblem::NavierStokesProblem(const InputParameters & parameters) : FEProblem(parameters)
47 : #if PETSC_RELEASE_LESS_THAN(3, 20, 0)
48 : {
49 : mooseError("The preconditioning techniques made available through this class require a PETSc "
50 : "version of at least 3.20");
51 : }
52 : #else
53 : ,
54 165 : _commute_lsc(getParam<bool>("commute_lsc")),
55 330 : _mass_matrix(getParam<TagName>("mass_matrix")),
56 330 : _L_matrix(getParam<TagName>("L_matrix")),
57 165 : _have_mass_matrix(!_mass_matrix.empty()),
58 165 : _have_L_matrix(!_L_matrix.empty()),
59 330 : _pressure_mass_matrix_as_pre(getParam<bool>("use_pressure_mass_matrix")),
60 495 : _schur_fs_index(getParam<std::vector<unsigned int>>("schur_fs_index"))
61 : {
62 165 : if (_commute_lsc)
63 : {
64 19 : if (!_have_mass_matrix)
65 0 : paramError("mass_matrix",
66 : "A pressure mass matrix must be provided if we are commuting the LSC commutator.");
67 19 : if (!_have_L_matrix)
68 0 : paramError("L_matrix",
69 : "A matrix corresponding to the viscous component of the momentum equation must be "
70 : "provided if we are commuting the LSC commutator.");
71 : }
72 146 : else if (_have_L_matrix)
73 0 : paramError("L_matrix",
74 : "If not commuting the LSC commutator, then the 'L_matrix' should not be provided "
75 : "because it will not be used. For Elman LSC preconditioning, L will be computed "
76 : "automatically using system matrix data (e.g. the off-diagonal blocks in the "
77 : "velocity-pressure system).");
78 :
79 165 : if (_pressure_mass_matrix_as_pre && !_have_mass_matrix)
80 0 : paramError("mass_matrix",
81 : "If 'use_pressure_mass_matrix', then a pressure 'mass_matrix' must be provided");
82 165 : }
83 :
84 330 : NavierStokesProblem::~NavierStokesProblem()
85 : {
86 : auto ierr = (PetscErrorCode)0;
87 165 : if (_Q_scale)
88 : {
89 116 : ierr = MatDestroy(&_Q_scale);
90 116 : CHKERRABORT(this->comm().get(), ierr);
91 : }
92 165 : if (_L)
93 : {
94 14 : ierr = MatDestroy(&_L);
95 14 : CHKERRABORT(this->comm().get(), ierr);
96 : }
97 330 : }
98 :
99 : KSP
100 3338 : NavierStokesProblem::findSchurKSP(KSP node, const unsigned int tree_position)
101 : {
102 3338 : auto it = _schur_fs_index.begin() + tree_position;
103 3338 : if (it == _schur_fs_index.end())
104 : return node;
105 :
106 : PC fs_pc;
107 : PetscInt num_splits;
108 : KSP * subksp;
109 : KSP next_ksp;
110 : IS is;
111 : PetscBool is_fs;
112 :
113 1432 : auto sub_ksp_index = *it;
114 :
115 : // Get the preconditioner associated with the linear solver
116 1432 : LibmeshPetscCall(KSPGetPC(node, &fs_pc));
117 :
118 : // Verify the preconditioner is a field split preconditioner
119 1432 : LibmeshPetscCall(PetscObjectTypeCompare((PetscObject)fs_pc, PCFIELDSPLIT, &is_fs));
120 1432 : if (!is_fs)
121 0 : mooseError("Not a field split. Please check the 'schur_fs_index' parameter");
122 :
123 : // Setup the preconditioner. We need to call this first in order to be able to retrieve the sub
124 : // ksps and sub index sets associated with the splits
125 1432 : LibmeshPetscCall(PCSetUp(fs_pc));
126 :
127 : // Get the linear solvers associated with each split
128 1432 : LibmeshPetscCall(PCFieldSplitGetSubKSP(fs_pc, &num_splits, &subksp));
129 1432 : next_ksp = subksp[sub_ksp_index];
130 :
131 : // Get the index set for the split at this level of the tree we are traversing to the Schur
132 : // complement preconditioner
133 1432 : LibmeshPetscCall(PCFieldSplitGetISByIndex(fs_pc, sub_ksp_index, &is));
134 :
135 : // Store this tree level's index set, which we will eventually use to get the sub-matrices
136 : // required for our preconditioning process from the system matrix
137 1432 : _index_sets.push_back(is);
138 :
139 : // Free the array of sub linear solvers that got allocated in the PCFieldSplitGetSubKSP call
140 1432 : LibmeshPetscCall(PetscFree(subksp));
141 :
142 : // Continue traversing down the tree towards the Schur complement linear solver/preconditioner
143 1432 : return findSchurKSP(next_ksp, tree_position + 1);
144 : }
145 :
146 : void
147 1906 : NavierStokesProblem::setupLSCMatrices(KSP schur_ksp)
148 : {
149 : KSP * subksp; // This will be length two, with the former being the A KSP and the latter being the
150 : // Schur complement KSP
151 : KSP schur_complement_ksp;
152 : PC schur_pc, lsc_pc;
153 : PetscInt num_splits;
154 : Mat lsc_pc_pmat;
155 : IS velocity_is, pressure_is;
156 : PetscInt rstart, rend;
157 : PetscBool is_lsc, is_fs;
158 : std::vector<Mat> intermediate_Qs;
159 : std::vector<Mat> intermediate_Ls;
160 :
161 : // Get the preconditioner for the linear solver. It must be a field split preconditioner
162 1906 : LibmeshPetscCall(KSPGetPC(schur_ksp, &schur_pc));
163 1906 : LibmeshPetscCall(PetscObjectTypeCompare((PetscObject)schur_pc, PCFIELDSPLIT, &is_fs));
164 1906 : if (!is_fs)
165 0 : mooseError("Not a field split. Please check the 'schur_fs_index' parameter");
166 :
167 : // The mass matrix. This will correspond to velocity degrees of freedom for Elman LSC and pressure
168 : // degrees of freedom for Olshanskii LSC or when directly using the mass matrix as a
169 : // preconditioner for the Schur complement
170 : Mat global_Q = nullptr;
171 1906 : if (_have_mass_matrix)
172 : global_Q =
173 1906 : static_cast<PetscMatrix<Number> &>(getNonlinearSystemBase(0).getMatrix(massMatrixTagID()))
174 1906 : .mat();
175 : // The Poisson operator matrix corresponding to the velocity degrees of freedom. This is only used
176 : // and is required for Olshanskii LSC preconditioning
177 : Mat global_L = nullptr;
178 1906 : if (_have_L_matrix)
179 : global_L =
180 42 : static_cast<PetscMatrix<Number> &>(getNonlinearSystemBase(0).getMatrix(LMatrixTagID()))
181 42 : .mat();
182 :
183 : //
184 : // Process down from our system matrix to the sub-matrix containing the velocity-pressure dofs for
185 : // which we are going to be doing the Schur complement preconditioning
186 : //
187 :
188 1948 : auto process_intermediate_mats = [this](auto & intermediate_mats, auto parent_mat)
189 : {
190 : mooseAssert(parent_mat, "This should be non-null");
191 1948 : intermediate_mats.resize(_index_sets.size());
192 3380 : for (const auto i : index_range(_index_sets))
193 : {
194 1432 : auto intermediate_is = _index_sets[i];
195 : Mat intermediate_mat;
196 1432 : LibmeshPetscCall(MatCreateSubMatrix(i == 0 ? parent_mat : intermediate_mats[i - 1],
197 : intermediate_is,
198 : intermediate_is,
199 : MAT_INITIAL_MATRIX,
200 : &intermediate_mat));
201 1432 : intermediate_mats[i] = intermediate_mat;
202 : }
203 1948 : return _index_sets.empty() ? parent_mat : intermediate_mats.back();
204 1906 : };
205 :
206 : Mat our_parent_Q = nullptr;
207 1906 : if (_have_mass_matrix)
208 1906 : our_parent_Q = process_intermediate_mats(intermediate_Qs, global_Q);
209 : Mat our_parent_L = nullptr;
210 1906 : if (_have_L_matrix)
211 42 : our_parent_L = process_intermediate_mats(intermediate_Ls, global_L);
212 :
213 : // Setup the preconditioner. We need to call this first in order to be able to retrieve the sub
214 : // ksps and sub index sets associated with the splits
215 1906 : LibmeshPetscCall(PCSetUp(schur_pc));
216 :
217 : // There are always two splits in a Schur complement split. The zeroth split is the split with the
218 : // on-diagonals, e.g. the velocity dofs. Here we retrive the velocity dofs/index set
219 1906 : LibmeshPetscCall(PCFieldSplitGetISByIndex(schur_pc, 0, &velocity_is));
220 :
221 : // Get the rows of the parent velocity-pressure matrix that our process owns
222 1906 : LibmeshPetscCall(MatGetOwnershipRange(our_parent_Q, &rstart, &rend));
223 :
224 1906 : if (_commute_lsc)
225 : {
226 : // If we're commuting LSC, e.g. doing Olshanskii, the user must have provided a Poisson operator
227 : // matrix
228 : mooseAssert(our_parent_L, "This should be non-null");
229 :
230 42 : if (!_L)
231 : // If this is our first time in this routine, then we create the matrix
232 14 : LibmeshPetscCall(
233 : MatCreateSubMatrix(our_parent_L, velocity_is, velocity_is, MAT_INITIAL_MATRIX, &_L));
234 : else
235 : // Else we reuse the matrix
236 28 : LibmeshPetscCall(
237 : MatCreateSubMatrix(our_parent_L, velocity_is, velocity_is, MAT_REUSE_MATRIX, &_L));
238 : }
239 :
240 : // Get the local index set complement corresponding to the pressure dofs from the velocity dofs
241 1906 : LibmeshPetscCall(ISComplement(velocity_is, rstart, rend, &pressure_is));
242 :
243 : auto create_q_scale_submat =
244 1906 : [our_parent_Q, this, velocity_is, pressure_is](const auto & mat_initialization)
245 : {
246 1906 : if (_commute_lsc || _pressure_mass_matrix_as_pre)
247 : {
248 : // If we are doing Olshanskii or we are using the pressure matrix directly as the
249 : // preconditioner (no LSC), then we must have access to a pressure mass matrix
250 : mooseAssert(our_parent_Q, "This should be non-null");
251 : // Create a sub-matrix corresponding to the pressure index set
252 292 : LibmeshPetscCall(MatCreateSubMatrix(
253 : our_parent_Q, pressure_is, pressure_is, mat_initialization, &_Q_scale));
254 : }
255 1614 : else if (_have_mass_matrix) // If we don't have a mass matrix and the user has requested scaling
256 : // then the diagonal of A will be used
257 : {
258 : // The user passed us a mass matrix tag; we better have been able to obtain a parent Q in that
259 : // case
260 : mooseAssert(our_parent_Q, "This should be non-null");
261 : // We are not commuting LSC, so we are doing Elman, and the user has passed us a mass matrix
262 : // tag. In this case we are creating a velocity mass matrix, so we use the velocity index set
263 1614 : LibmeshPetscCall(MatCreateSubMatrix(
264 : our_parent_Q, velocity_is, velocity_is, mat_initialization, &_Q_scale));
265 : }
266 1906 : };
267 :
268 1906 : if (!_Q_scale)
269 : // We haven't allocated the scaling matrix yet
270 116 : create_q_scale_submat(MAT_INITIAL_MATRIX);
271 : else
272 : // We have allocated the scaling matrix, so we can reuse
273 1790 : create_q_scale_submat(MAT_REUSE_MATRIX);
274 :
275 : // We don't need the pressure index set anymore
276 1906 : LibmeshPetscCall(ISDestroy(&pressure_is));
277 :
278 : // Nor the intermediate matrices
279 3338 : for (auto & mat : intermediate_Qs)
280 1432 : LibmeshPetscCall(MatDestroy(&mat));
281 1906 : for (auto & mat : intermediate_Ls)
282 0 : LibmeshPetscCall(MatDestroy(&mat));
283 :
284 : // Get the sub KSP for the Schur split that corresponds to the linear solver for the Schur
285 : // complement (e.g. rank equivalent to the pressure rank)
286 1906 : LibmeshPetscCall(PCFieldSplitGetSubKSP(schur_pc, &num_splits, &subksp));
287 1906 : if (num_splits != 2)
288 0 : mooseError("The number of splits should be two");
289 : // The Schur complement linear solver is always at the first index (for the pressure dofs;
290 : // velocity dof KSP is at index 0)
291 1906 : schur_complement_ksp = subksp[1];
292 :
293 1906 : if (_pressure_mass_matrix_as_pre)
294 : {
295 : mooseAssert(_Q_scale, "This should be non-null");
296 : Mat S;
297 : // Set the Schur complement preconditioner to be the pressure mass matrix
298 250 : LibmeshPetscCall(PCFieldSplitSetSchurPre(schur_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, _Q_scale));
299 : // Get the Schur complement operator S, which in generic KSP speak is used for the operator A
300 250 : LibmeshPetscCall(KSPGetOperators(schur_complement_ksp, &S, NULL));
301 : // Set, in generic KSP speak, the operators A and P respectively. So our pressure mass matrix is
302 : // P
303 250 : LibmeshPetscCall(KSPSetOperators(schur_complement_ksp, S, _Q_scale));
304 : }
305 : else // We are doing LSC preconditioning
306 : {
307 : // Get the least squares commutator preconditioner for the Schur complement
308 1656 : LibmeshPetscCall(KSPGetPC(schur_complement_ksp, &lsc_pc));
309 : // Verify that it's indeed an LSC preconditioner
310 1656 : LibmeshPetscCall(PetscObjectTypeCompare(PetscObject(lsc_pc), PCLSC, &is_lsc));
311 1656 : if (!is_lsc)
312 0 : mooseError("Not an LSC PC. Please check the 'schur_fs_index' parameter");
313 :
314 : // Get the LSC preconditioner
315 1656 : LibmeshPetscCall(PCGetOperators(lsc_pc, NULL, &lsc_pc_pmat));
316 :
317 1656 : if (_commute_lsc)
318 : {
319 : // We're doing Olshanskii. We must have a user-provided Poisson operator
320 : mooseAssert(_L, "This should be non-null");
321 : // Attach our L matrix to the PETSc object. PETSc will use this during the preconditioner
322 : // application
323 42 : LibmeshPetscCall(PetscObjectCompose((PetscObject)lsc_pc_pmat, "LSC_L", (PetscObject)_L));
324 : // Olshanskii preconditioning requires a pressure mass matrix
325 : mooseAssert(_have_mass_matrix, "This is to verify we will enter the next conditional");
326 : }
327 1656 : if (_have_mass_matrix)
328 : {
329 : mooseAssert(_Q_scale, "This should be non-null");
330 : // Attach our scaling/mass matrix to the PETSc object. PETSc will use this during the
331 : // preconditioner application
332 1656 : LibmeshPetscCall(
333 : PetscObjectCompose((PetscObject)lsc_pc_pmat, "LSC_Qscale", (PetscObject)_Q_scale));
334 : }
335 : }
336 :
337 : // Free the sub-KSP array that was allocated during PCFieldSplitGetSubKSP
338 1906 : LibmeshPetscCall(PetscFree(subksp));
339 1906 : }
340 :
341 : PetscErrorCode
342 1906 : navierStokesKSPPreSolve(KSP root_ksp, Vec /*rhs*/, Vec /*x*/, void * context)
343 : {
344 : PetscFunctionBegin;
345 :
346 : auto * ns_problem = static_cast<NavierStokesProblem *>(context);
347 : ns_problem->clearIndexSets();
348 1906 : auto schur_ksp = ns_problem->findSchurKSP(root_ksp, 0);
349 1906 : ns_problem->setupLSCMatrices(schur_ksp);
350 :
351 1906 : PetscFunctionReturn(PETSC_SUCCESS);
352 : }
353 :
354 : void
355 370 : NavierStokesProblem::initPetscOutputAndSomeSolverSettings()
356 : {
357 370 : FEProblem::initPetscOutputAndSomeSolverSettings();
358 :
359 370 : if (!_have_mass_matrix)
360 : {
361 : mooseAssert(
362 : !_have_L_matrix,
363 : "If we don't have a mass matrix, which is only supported for traditional LSC "
364 : "preconditioning (e.g. Elman, not Olshanskii), then we also shouldn't have an L matrix "
365 : "because we automatically form the L matrix when doing traditional LSC preconditioning");
366 14 : return;
367 : }
368 :
369 : // Set the pre-KSP solve callback. At that time we will setup our Schur complement preconditioning
370 : KSP ksp;
371 356 : auto snes = currentNonlinearSystem().getSNES();
372 356 : LibmeshPetscCall(SNESGetKSP(snes, &ksp));
373 356 : LibmeshPetscCall(KSPSetPreSolve(ksp, &navierStokesKSPPreSolve, this));
374 : }
375 :
376 : #endif
|
