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