Actual source code: ex4.c
1: static char help[] = "Applies the 2023 preconditioner of Benzi and Faccio\n\n";
3: #include <petscmat.h>
4: #include <petscvec.h>
5: #include <petscis.h>
6: #include <petscksp.h>
8: /*
9: * This example reproduces the preconditioner outlined in Benzi's paper
10: * https://doi.org/10.1137/22M1505529. The problem considered is:
11: *
12: * (A + gamma UU^T)x = b
13: *
14: * whose structure arises from, for example, grad-div stabilization in the
15: * Navier-Stokes momentum equation. In the code we will also refer to
16: * gamma UU^T as J. The preconditioner developed by Benzi is:
17: *
18: * P_alpha = (A + alpha I)(alpha I + gamma UU^T)
19: *
20: * Another variant which may yield better convergence depending on the specific
21: * problem is
22: *
23: * P_alpha = (A + alpha D) D^-1 (alpha D + gamma UU^T)
24: *
25: * where D = diag(A + gamma UU^T). This is the variant implemented
26: * here. Application of the preconditioner involves (approximate) solution of
27: * two systems, one with (A + alpha D), and another with (alpha D + gamma
28: * UU^T). For small alpha (which generally yields the best overall
29: * preconditioner), (alpha D + gamma UU^T) is ill-conditioned. To combat this we
30: * solve (alpha D + gamma UU^T) using the Sherman-Morrison-Woodbury (SMW) matrix
31: * identity, which effectively converts the grad-div structure to a much nicer
32: * div-grad (laplacian) structure.
33: *
34: * The matrices used as input can be generated by running the matlab/octave
35: * program IFISS. The particular matrices checked into the datafiles repository
36: * and used in testing of this example correspond to a leaky lid-driven cavity
37: * with a stretched grid and Q2-Q1 finite elements. The matrices are taken from
38: * the last iteration of a Picard solve with tolerance 1e-8 with a viscosity of
39: * 0.1 and a 32x32 grid. We summarize below iteration counts from running this
40: * preconditioner for different grids and viscosity with a KSP tolerance of 1e-6.
41: *
42: * 32x32 64x64 128x128
43: * 0.1 28 36 43
44: * 0.01 59 75 73
45: * 0.002 136 161 167
46: *
47: * A reader of Benzi's paper will note that the performance shown above with
48: * respect to decreasing viscosity is significantly worse than in the
49: * paper. This is actually because of the choice of RHS. In Benzi's work, the
50: * RHS was generated by multiplying the operator with a vector of 1s whereas
51: * here we generate the RHS using a random vector. The iteration counts from the
52: * Benzi paper can be reproduced by changing the RHS generation in this example,
53: * but we choose to use the more difficult RHS as the resulting performance may
54: * more closely match what users experience in "physical" contexts.
55: */
57: PetscErrorCode CreateAndLoadMat(const char *mat_name, Mat *mat)
58: {
59: PetscViewer viewer;
60: char file[PETSC_MAX_PATH_LEN];
61: char flag_name[10] = "-f";
62: PetscBool flg;
64: PetscFunctionBeginUser;
65: PetscCall(PetscOptionsGetString(NULL, NULL, strcat(flag_name, mat_name), file, sizeof(file), &flg));
66: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate file with the -f<mat_name> option");
67: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &viewer));
68: PetscCall(MatCreate(PETSC_COMM_WORLD, mat));
69: PetscCall(MatSetType(*mat, MATAIJ));
70: PetscCall(PetscObjectSetName((PetscObject)*mat, mat_name));
71: PetscCall(MatSetFromOptions(*mat));
72: PetscCall(MatLoad(*mat, viewer));
73: PetscCall(PetscViewerDestroy(&viewer));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: typedef struct {
78: Mat U, UT, D, aD, aDinv, I_plus_gammaUTaDinvU;
79: PC smw_cholesky;
80: PetscReal gamma, alpha;
81: PetscBool setup_called;
82: } SmwPCCtx;
84: PetscErrorCode SmwSetup(PC pc)
85: {
86: SmwPCCtx *ctx;
87: Vec aDVec;
89: PetscFunctionBeginUser;
90: PetscCall(PCShellGetContext(pc, &ctx));
92: if (ctx->setup_called) PetscFunctionReturn(PETSC_SUCCESS);
94: // Create aD
95: PetscCall(MatDuplicate(ctx->D, MAT_COPY_VALUES, &ctx->aD));
96: PetscCall(MatScale(ctx->aD, ctx->alpha));
98: // Create aDinv
99: PetscCall(MatDuplicate(ctx->aD, MAT_DO_NOT_COPY_VALUES, &ctx->aDinv));
100: PetscCall(MatCreateVecs(ctx->aD, &aDVec, NULL));
101: PetscCall(MatGetDiagonal(ctx->aD, aDVec));
102: PetscCall(VecReciprocal(aDVec));
103: PetscCall(MatDiagonalSet(ctx->aDinv, aDVec, INSERT_VALUES));
105: // Create UT
106: PetscCall(MatTranspose(ctx->U, MAT_INITIAL_MATRIX, &ctx->UT));
108: // Create sum Mat
109: PetscCall(MatMatMatMult(ctx->UT, ctx->aDinv, ctx->U, MAT_INITIAL_MATRIX, PETSC_CURRENT, &ctx->I_plus_gammaUTaDinvU));
110: PetscCall(MatScale(ctx->I_plus_gammaUTaDinvU, ctx->gamma));
111: PetscCall(MatShift(ctx->I_plus_gammaUTaDinvU, 1.));
113: PetscCall(PCCreate(PETSC_COMM_WORLD, &ctx->smw_cholesky));
114: PetscCall(PCSetType(ctx->smw_cholesky, PCCHOLESKY));
115: PetscCall(PCSetOperators(ctx->smw_cholesky, ctx->I_plus_gammaUTaDinvU, ctx->I_plus_gammaUTaDinvU));
116: PetscCall(PCSetOptionsPrefix(ctx->smw_cholesky, "smw_"));
117: PetscCall(PCSetFromOptions(ctx->smw_cholesky));
118: PetscCall(PCSetUp(ctx->smw_cholesky));
120: PetscCall(VecDestroy(&aDVec));
122: ctx->setup_called = PETSC_TRUE;
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: PetscErrorCode SmwApply(PC pc, Vec x, Vec y)
127: {
128: SmwPCCtx *ctx;
129: Vec vel0, pressure0, pressure1;
131: PetscFunctionBeginUser;
132: PetscCall(PCShellGetContext(pc, &ctx));
134: PetscCall(MatCreateVecs(ctx->UT, &vel0, &pressure0));
135: PetscCall(VecDuplicate(pressure0, &pressure1));
137: // First term
138: PetscCall(MatMult(ctx->aDinv, x, vel0));
139: PetscCall(MatMult(ctx->UT, vel0, pressure0));
140: PetscCall(PCApply(ctx->smw_cholesky, pressure0, pressure1));
141: PetscCall(MatMult(ctx->U, pressure1, vel0));
142: PetscCall(MatMult(ctx->aDinv, vel0, y));
143: PetscCall(VecScale(y, -ctx->gamma));
145: // Second term
146: PetscCall(MatMult(ctx->aDinv, x, vel0));
148: PetscCall(VecAXPY(y, 1, vel0));
150: PetscCall(VecDestroy(&vel0));
151: PetscCall(VecDestroy(&pressure0));
152: PetscCall(VecDestroy(&pressure1));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: int main(int argc, char **args)
157: {
158: Mat A, B, Q, Acondensed, Bcondensed, BT, J, AplusJ, QInv, D, AplusD, JplusD, U;
159: Mat AplusJarray[2];
160: Vec bound, x, b, Qdiag, DVec;
161: PetscBool flg;
162: PetscViewer viewer;
163: char file[PETSC_MAX_PATH_LEN];
164: PetscInt *boundary_indices;
165: PetscInt boundary_indices_size, am, an, bm, bn, condensed_am, astart, aend, Dstart, Dend, num_local_bnd_dofs = 0;
166: const PetscScalar zero = 0;
167: IS boundary_is, bulk_is;
168: KSP ksp;
169: PC pc, pcA, pcJ;
170: PetscRandom rctx;
171: PetscReal *boundary_indices_values;
172: PetscReal gamma = 100, alpha = .01;
173: PetscMPIInt rank;
174: SmwPCCtx ctx;
176: PetscFunctionBeginUser;
177: PetscCall(PetscInitialize(&argc, &args, NULL, help));
178: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
180: PetscCall(CreateAndLoadMat("A", &A));
181: PetscCall(CreateAndLoadMat("B", &B));
182: PetscCall(CreateAndLoadMat("Q", &Q));
184: PetscCall(PetscOptionsGetString(NULL, NULL, "-fbound", file, sizeof(file), &flg));
185: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate file with the -fbound option");
187: if (rank == 0) {
188: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, file, FILE_MODE_READ, &viewer));
189: PetscCall(VecCreate(PETSC_COMM_SELF, &bound));
190: PetscCall(PetscObjectSetName((PetscObject)bound, "bound"));
191: PetscCall(VecSetType(bound, VECSEQ));
192: PetscCall(VecLoad(bound, viewer));
193: PetscCall(PetscViewerDestroy(&viewer));
194: PetscCall(VecGetLocalSize(bound, &boundary_indices_size));
195: PetscCall(VecGetArray(bound, &boundary_indices_values));
196: }
197: PetscCallMPI(MPI_Bcast(&boundary_indices_size, 1, MPIU_INT, 0, PETSC_COMM_WORLD));
198: if (rank != 0) PetscCall(PetscMalloc1(boundary_indices_size, &boundary_indices_values));
199: PetscCallMPI(MPI_Bcast(boundary_indices_values, (PetscMPIInt)boundary_indices_size, MPIU_SCALAR, 0, PETSC_COMM_WORLD));
201: PetscCall(MatGetSize(A, &am, NULL));
202: // The total number of dofs for a given velocity component
203: const PetscInt nc = am / 2;
204: PetscCall(MatGetOwnershipRange(A, &astart, &aend));
206: PetscCall(PetscMalloc1(2 * boundary_indices_size, &boundary_indices));
208: //
209: // The dof index ordering appears to be all vx dofs and then all vy dofs.
210: //
212: // First do vx
213: for (PetscInt i = 0; i < boundary_indices_size; ++i) {
214: // MATLAB uses 1-based indexing
215: const PetscInt bnd_dof = (PetscInt)boundary_indices_values[i] - 1;
216: if (bnd_dof >= astart && bnd_dof < aend) boundary_indices[num_local_bnd_dofs++] = bnd_dof;
217: }
219: // Now vy
220: for (PetscInt i = 0; i < boundary_indices_size; ++i) {
221: // MATLAB uses 1-based indexing
222: const PetscInt bnd_dof = ((PetscInt)boundary_indices_values[i] - 1) + nc;
223: if (bnd_dof >= astart && bnd_dof < aend) boundary_indices[num_local_bnd_dofs++] = bnd_dof;
224: }
225: if (rank == 0) PetscCall(VecRestoreArray(bound, &boundary_indices_values));
226: else PetscCall(PetscFree(boundary_indices_values));
228: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, num_local_bnd_dofs, boundary_indices, PETSC_USE_POINTER, &boundary_is));
229: PetscCall(ISComplement(boundary_is, astart, aend, &bulk_is));
231: PetscCall(MatCreateSubMatrix(A, bulk_is, bulk_is, MAT_INITIAL_MATRIX, &Acondensed));
232: // Can't pass null for row index set :-(
233: PetscCall(MatTranspose(B, MAT_INPLACE_MATRIX, &B));
234: PetscCall(MatCreateSubMatrix(B, bulk_is, NULL, MAT_INITIAL_MATRIX, &Bcondensed));
235: PetscCall(MatGetLocalSize(Acondensed, &am, &an));
236: PetscCall(MatGetLocalSize(Bcondensed, &bm, &bn));
238: // Create QInv
239: PetscCall(MatCreateVecs(Q, &Qdiag, NULL));
240: PetscCall(MatGetDiagonal(Q, Qdiag));
241: PetscCall(VecReciprocal(Qdiag));
242: PetscCall(MatDuplicate(Q, MAT_DO_NOT_COPY_VALUES, &QInv));
243: PetscCall(MatDiagonalSet(QInv, Qdiag, INSERT_VALUES));
244: PetscCall(MatAssemblyBegin(QInv, MAT_FINAL_ASSEMBLY));
245: PetscCall(MatAssemblyEnd(QInv, MAT_FINAL_ASSEMBLY));
247: // Create J
248: PetscCall(MatTranspose(Bcondensed, MAT_INITIAL_MATRIX, &BT));
249: PetscCall(MatMatMatMult(Bcondensed, QInv, BT, MAT_INITIAL_MATRIX, PETSC_CURRENT, &J));
250: PetscCall(MatScale(J, gamma));
252: // Create sum of A + J
253: AplusJarray[0] = Acondensed;
254: AplusJarray[1] = J;
255: PetscCall(MatCreateComposite(PETSC_COMM_WORLD, 2, AplusJarray, &AplusJ));
257: // Create decomposition matrices
258: // We've already used Qdiag, which currently represents Q^-1, for its necessary purposes. Let's
259: // convert it to represent Q^(-1/2)
260: PetscCall(VecSqrtAbs(Qdiag));
261: // We can similarly reuse Qinv
262: PetscCall(MatDiagonalSet(QInv, Qdiag, INSERT_VALUES));
263: PetscCall(MatAssemblyBegin(QInv, MAT_FINAL_ASSEMBLY));
264: PetscCall(MatAssemblyEnd(QInv, MAT_FINAL_ASSEMBLY));
265: // Create U
266: PetscCall(MatMatMult(Bcondensed, QInv, MAT_INITIAL_MATRIX, PETSC_CURRENT, &U));
268: // Create x and b
269: PetscCall(MatCreateVecs(AplusJ, &x, &b));
270: PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx));
271: PetscCall(VecSetRandom(x, rctx));
272: PetscCall(PetscRandomDestroy(&rctx));
273: PetscCall(MatMult(AplusJ, x, b));
275: // Compute preconditioner operators
276: PetscCall(MatGetLocalSize(Acondensed, &condensed_am, NULL));
277: PetscCall(MatCreate(PETSC_COMM_WORLD, &D));
278: PetscCall(MatSetType(D, MATAIJ));
279: PetscCall(MatSetSizes(D, condensed_am, condensed_am, PETSC_DETERMINE, PETSC_DETERMINE));
280: PetscCall(MatGetOwnershipRange(D, &Dstart, &Dend));
281: for (PetscInt i = Dstart; i < Dend; ++i) PetscCall(MatSetValues(D, 1, &i, 1, &i, &zero, INSERT_VALUES));
282: PetscCall(MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY));
283: PetscCall(MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY));
284: PetscCall(MatCreateVecs(D, &DVec, NULL));
285: PetscCall(MatGetDiagonal(AplusJ, DVec));
286: PetscCall(MatDiagonalSet(D, DVec, INSERT_VALUES));
287: PetscCall(MatDuplicate(Acondensed, MAT_COPY_VALUES, &AplusD));
288: PetscCall(MatAXPY(AplusD, alpha, D, SUBSET_NONZERO_PATTERN));
289: PetscCall(MatDuplicate(J, MAT_COPY_VALUES, &JplusD));
290: PetscCall(MatAXPY(JplusD, alpha, D, SUBSET_NONZERO_PATTERN));
292: // Initialize our SMW context
293: ctx.U = U;
294: ctx.D = D;
295: ctx.gamma = gamma;
296: ctx.alpha = alpha;
297: ctx.setup_called = PETSC_FALSE;
299: // Set preconditioner operators
300: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
301: PetscCall(KSPSetType(ksp, KSPFGMRES));
302: PetscCall(KSPSetOperators(ksp, AplusJ, AplusJ));
303: PetscCall(KSPSetNormType(ksp, KSP_NORM_UNPRECONDITIONED));
304: PetscCall(KSPGMRESSetRestart(ksp, 300));
305: PetscCall(KSPGetPC(ksp, &pc));
306: PetscCall(PCSetType(pc, PCCOMPOSITE));
307: PetscCall(PCCompositeSetType(pc, PC_COMPOSITE_SPECIAL));
308: PetscCall(PCCompositeAddPCType(pc, PCILU));
309: PetscCall(PCCompositeAddPCType(pc, PCSHELL));
310: PetscCall(PCCompositeGetPC(pc, 0, &pcA));
311: PetscCall(PCCompositeGetPC(pc, 1, &pcJ));
312: PetscCall(PCSetOperators(pcA, AplusD, AplusD));
313: PetscCall(PCSetOperators(pcJ, JplusD, JplusD));
314: PetscCall(PCShellSetContext(pcJ, &ctx));
315: PetscCall(PCShellSetApply(pcJ, SmwApply));
316: PetscCall(PCShellSetSetUp(pcJ, SmwSetup));
317: PetscCall(PCCompositeSpecialSetAlphaMat(pc, D));
319: // Solve
320: PetscCall(KSPSetFromOptions(ksp));
321: PetscCall(KSPSolve(ksp, b, x));
323: PetscCall(MatDestroy(&A));
324: PetscCall(MatDestroy(&B));
325: PetscCall(MatDestroy(&Q));
326: PetscCall(MatDestroy(&Acondensed));
327: PetscCall(MatDestroy(&Bcondensed));
328: PetscCall(MatDestroy(&BT));
329: PetscCall(MatDestroy(&J));
330: PetscCall(MatDestroy(&AplusJ));
331: PetscCall(MatDestroy(&QInv));
332: PetscCall(MatDestroy(&D));
333: PetscCall(MatDestroy(&AplusD));
334: PetscCall(MatDestroy(&JplusD));
335: PetscCall(MatDestroy(&U));
336: if (rank == 0) PetscCall(VecDestroy(&bound));
337: PetscCall(VecDestroy(&x));
338: PetscCall(VecDestroy(&b));
339: PetscCall(VecDestroy(&Qdiag));
340: PetscCall(VecDestroy(&DVec));
341: PetscCall(ISDestroy(&boundary_is));
342: PetscCall(ISDestroy(&bulk_is));
343: PetscCall(KSPDestroy(&ksp));
344: PetscCall(PetscFree(boundary_indices));
345: PetscCall(MatDestroy(&ctx.UT));
346: PetscCall(MatDestroy(&ctx.I_plus_gammaUTaDinvU));
347: PetscCall(MatDestroy(&ctx.aD));
348: PetscCall(MatDestroy(&ctx.aDinv));
349: PetscCall(PCDestroy(&ctx.smw_cholesky));
351: PetscCall(PetscFinalize());
352: return 0;
353: }
355: /*TEST
357: build:
358: requires: !complex
360: test:
361: args: -fA ${DATAFILESPATH}/matrices/ifiss/A -fB ${DATAFILESPATH}/matrices/ifiss/B -fQ ${DATAFILESPATH}/matrices/ifiss/Q -fbound ${DATAFILESPATH}/is/ifiss/bound -ksp_monitor
362: requires: datafilespath defined(PETSC_USE_64BIT_INDICES) !complex double
364: test:
365: suffix: 2
366: nsize: 2
367: args: -fA ${DATAFILESPATH}/matrices/ifiss/A -fB ${DATAFILESPATH}/matrices/ifiss/B -fQ ${DATAFILESPATH}/matrices/ifiss/Q -fbound ${DATAFILESPATH}/is/ifiss/bound -ksp_monitor
368: requires: datafilespath defined(PETSC_USE_64BIT_INDICES) !complex double strumpack
370: TEST*/