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