Actual source code: ex11.c
1: static const char help[] = "Solves a Q1-P0 Stokes problem from Underworld.\n\
2: \n\
3: You can obtain a sample matrix from http://ftp.mcs.anl.gov/pub/petsc/Datafiles/matrices/underworld32.gz\n\
4: and run with -f underworld32.gz\n\n";
6: #include <petscksp.h>
7: #include <petscdmda.h>
9: static PetscErrorCode replace_submats(Mat A, IS isu, IS isp)
10: {
11: Mat A11, A22, A12, A21;
12: Mat nA11, nA22, nA12, nA21;
13: const char *prefix;
16: MatCreateSubMatrix(A, isu, isu, MAT_INITIAL_MATRIX, &A11);
17: MatCreateSubMatrix(A, isu, isp, MAT_INITIAL_MATRIX, &A12);
18: MatCreateSubMatrix(A, isp, isu, MAT_INITIAL_MATRIX, &A21);
19: MatCreateSubMatrix(A, isp, isp, MAT_INITIAL_MATRIX, &A22);
20: MatDuplicate(A11, MAT_COPY_VALUES, &nA11);
21: MatDuplicate(A12, MAT_COPY_VALUES, &nA12);
22: MatDuplicate(A21, MAT_COPY_VALUES, &nA21);
23: MatDuplicate(A22, MAT_COPY_VALUES, &nA22);
24: MatGetOptionsPrefix(A11, &prefix);
25: MatSetOptionsPrefix(nA11, prefix);
26: MatGetOptionsPrefix(A22, &prefix);
27: MatSetOptionsPrefix(nA22, prefix);
28: MatNestSetSubMat(A, 0, 0, nA11);
29: MatNestSetSubMat(A, 0, 1, nA12);
30: MatNestSetSubMat(A, 1, 0, nA21);
31: MatNestSetSubMat(A, 1, 1, nA22);
32: MatDestroy(&A11);
33: MatDestroy(&A12);
34: MatDestroy(&A21);
35: MatDestroy(&A22);
36: MatDestroy(&nA11);
37: MatDestroy(&nA12);
38: MatDestroy(&nA21);
39: MatDestroy(&nA22);
40: return 0;
41: }
43: PetscErrorCode LSCLoadTestOperators(Mat *A11, Mat *A12, Mat *A21, Mat *A22, Vec *b1, Vec *b2)
44: {
45: PetscViewer viewer;
46: char filename[PETSC_MAX_PATH_LEN];
47: PetscBool flg;
50: MatCreate(PETSC_COMM_WORLD, A11);
51: MatCreate(PETSC_COMM_WORLD, A12);
52: MatCreate(PETSC_COMM_WORLD, A21);
53: MatCreate(PETSC_COMM_WORLD, A22);
54: MatSetOptionsPrefix(*A11, "a11_");
55: MatSetOptionsPrefix(*A22, "a22_");
56: MatSetFromOptions(*A11);
57: MatSetFromOptions(*A22);
58: VecCreate(PETSC_COMM_WORLD, b1);
59: VecCreate(PETSC_COMM_WORLD, b2);
60: /* Load matrices from a Q1-P0 discretisation of variable viscosity Stokes. The matrix blocks are packed into one file. */
61: PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg);
63: PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer);
64: MatLoad(*A11, viewer);
65: MatLoad(*A12, viewer);
66: MatLoad(*A21, viewer);
67: MatLoad(*A22, viewer);
68: VecLoad(*b1, viewer);
69: VecLoad(*b2, viewer);
70: PetscViewerDestroy(&viewer);
71: return 0;
72: }
74: PetscErrorCode LoadTestMatrices(Mat *_A, Vec *_x, Vec *_b, IS *_isu, IS *_isp)
75: {
76: Vec f, h, x, b, bX[2];
77: Mat A, Auu, Aup, Apu, App, bA[2][2];
78: IS is_u, is_p, bis[2];
79: PetscInt lnu, lnp, nu, np, i, start_u, end_u, start_p, end_p;
80: VecScatter *vscat;
81: PetscMPIInt rank;
84: /* fetch test matrices and vectors */
85: LSCLoadTestOperators(&Auu, &Aup, &Apu, &App, &f, &h);
87: /* build the mat-nest */
88: VecGetSize(f, &nu);
89: VecGetSize(h, &np);
91: VecGetLocalSize(f, &lnu);
92: VecGetLocalSize(h, &lnp);
94: VecGetOwnershipRange(f, &start_u, &end_u);
95: VecGetOwnershipRange(h, &start_p, &end_p);
97: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
98: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] lnu = %" PetscInt_FMT " | lnp = %" PetscInt_FMT " \n", rank, lnu, lnp);
99: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] s_u = %" PetscInt_FMT " | e_u = %" PetscInt_FMT " \n", rank, start_u, end_u);
100: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] s_p = %" PetscInt_FMT " | e_p = %" PetscInt_FMT " \n", rank, start_p, end_p);
101: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] is_u (offset) = %" PetscInt_FMT " \n", rank, start_u + start_p);
102: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] is_p (offset) = %" PetscInt_FMT " \n", rank, start_u + start_p + lnu);
103: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
105: ISCreateStride(PETSC_COMM_WORLD, lnu, start_u + start_p, 1, &is_u);
106: ISCreateStride(PETSC_COMM_WORLD, lnp, start_u + start_p + lnu, 1, &is_p);
108: bis[0] = is_u;
109: bis[1] = is_p;
110: bA[0][0] = Auu;
111: bA[0][1] = Aup;
112: bA[1][0] = Apu;
113: bA[1][1] = App;
114: MatCreateNest(PETSC_COMM_WORLD, 2, bis, 2, bis, &bA[0][0], &A);
115: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
116: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
118: /* Pull f,h into b */
119: MatCreateVecs(A, &b, &x);
120: bX[0] = f;
121: bX[1] = h;
122: PetscMalloc1(2, &vscat);
123: for (i = 0; i < 2; i++) {
124: VecScatterCreate(b, bis[i], bX[i], NULL, &vscat[i]);
125: VecScatterBegin(vscat[i], bX[i], b, INSERT_VALUES, SCATTER_REVERSE);
126: VecScatterEnd(vscat[i], bX[i], b, INSERT_VALUES, SCATTER_REVERSE);
127: VecScatterDestroy(&vscat[i]);
128: }
130: PetscFree(vscat);
131: MatDestroy(&Auu);
132: MatDestroy(&Aup);
133: MatDestroy(&Apu);
134: MatDestroy(&App);
135: VecDestroy(&f);
136: VecDestroy(&h);
138: *_isu = is_u;
139: *_isp = is_p;
140: *_A = A;
141: *_x = x;
142: *_b = b;
143: return 0;
144: }
146: PetscErrorCode port_lsd_bfbt(void)
147: {
148: Mat A, P;
149: Vec x, b;
150: KSP ksp_A;
151: PC pc_A;
152: IS isu, isp;
153: PetscBool test_fs = PETSC_FALSE;
156: LoadTestMatrices(&A, &x, &b, &isu, &isp);
157: PetscOptionsGetBool(NULL, NULL, "-test_fs", &test_fs, NULL);
158: if (!test_fs) {
159: PetscObjectReference((PetscObject)A);
160: P = A;
161: } else {
162: MatDuplicate(A, MAT_COPY_VALUES, &P);
163: }
164: KSPCreate(PETSC_COMM_WORLD, &ksp_A);
165: KSPSetOptionsPrefix(ksp_A, "fc_");
166: KSPSetOperators(ksp_A, A, P);
168: KSPSetFromOptions(ksp_A);
169: KSPGetPC(ksp_A, &pc_A);
170: PetscObjectTypeCompare((PetscObject)pc_A, PCLU, &test_fs);
171: if (test_fs) {
172: MatDestroy(&P);
173: MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &P);
174: KSPSetOperators(ksp_A, A, P);
175: KSPSetFromOptions(ksp_A);
176: KSPSolve(ksp_A, b, x);
177: } else {
178: PCFieldSplitSetBlockSize(pc_A, 2);
179: PCFieldSplitSetIS(pc_A, "velocity", isu);
180: PCFieldSplitSetIS(pc_A, "pressure", isp);
181: KSPSolve(ksp_A, b, x);
183: /* Pull u,p out of x */
184: {
185: PetscInt loc;
186: PetscReal max, norm;
187: PetscScalar sum;
188: Vec uvec, pvec;
189: VecScatter uscat, pscat;
190: Mat A11, A22;
192: /* grab matrices and create the compatible u,p vectors */
193: MatCreateSubMatrix(A, isu, isu, MAT_INITIAL_MATRIX, &A11);
194: MatCreateSubMatrix(A, isp, isp, MAT_INITIAL_MATRIX, &A22);
196: MatCreateVecs(A11, &uvec, NULL);
197: MatCreateVecs(A22, &pvec, NULL);
199: /* perform the scatter from x -> (u,p) */
200: VecScatterCreate(x, isu, uvec, NULL, &uscat);
201: VecScatterBegin(uscat, x, uvec, INSERT_VALUES, SCATTER_FORWARD);
202: VecScatterEnd(uscat, x, uvec, INSERT_VALUES, SCATTER_FORWARD);
204: VecScatterCreate(x, isp, pvec, NULL, &pscat);
205: VecScatterBegin(pscat, x, pvec, INSERT_VALUES, SCATTER_FORWARD);
206: VecScatterEnd(pscat, x, pvec, INSERT_VALUES, SCATTER_FORWARD);
208: PetscPrintf(PETSC_COMM_WORLD, "-- vector vector values --\n");
209: VecMin(uvec, &loc, &max);
210: PetscPrintf(PETSC_COMM_WORLD, " Min(u) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc);
211: VecMax(uvec, &loc, &max);
212: PetscPrintf(PETSC_COMM_WORLD, " Max(u) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc);
213: VecNorm(uvec, NORM_2, &norm);
214: PetscPrintf(PETSC_COMM_WORLD, " Norm(u) = %1.6f \n", (double)norm);
215: VecSum(uvec, &sum);
216: PetscPrintf(PETSC_COMM_WORLD, " Sum(u) = %1.6f \n", (double)PetscRealPart(sum));
218: PetscPrintf(PETSC_COMM_WORLD, "-- pressure vector values --\n");
219: VecMin(pvec, &loc, &max);
220: PetscPrintf(PETSC_COMM_WORLD, " Min(p) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc);
221: VecMax(pvec, &loc, &max);
222: PetscPrintf(PETSC_COMM_WORLD, " Max(p) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc);
223: VecNorm(pvec, NORM_2, &norm);
224: PetscPrintf(PETSC_COMM_WORLD, " Norm(p) = %1.6f \n", (double)norm);
225: VecSum(pvec, &sum);
226: PetscPrintf(PETSC_COMM_WORLD, " Sum(p) = %1.6f \n", (double)PetscRealPart(sum));
228: PetscPrintf(PETSC_COMM_WORLD, "-- Full vector values --\n");
229: VecMin(x, &loc, &max);
230: PetscPrintf(PETSC_COMM_WORLD, " Min(u,p) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc);
231: VecMax(x, &loc, &max);
232: PetscPrintf(PETSC_COMM_WORLD, " Max(u,p) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc);
233: VecNorm(x, NORM_2, &norm);
234: PetscPrintf(PETSC_COMM_WORLD, " Norm(u,p) = %1.6f \n", (double)norm);
235: VecSum(x, &sum);
236: PetscPrintf(PETSC_COMM_WORLD, " Sum(u,p) = %1.6f \n", (double)PetscRealPart(sum));
238: VecScatterDestroy(&uscat);
239: VecScatterDestroy(&pscat);
240: VecDestroy(&uvec);
241: VecDestroy(&pvec);
242: MatDestroy(&A11);
243: MatDestroy(&A22);
244: }
246: /* test second solve by changing the mat associated to the MATNEST blocks */
247: {
248: replace_submats(A, isu, isp);
249: replace_submats(P, isu, isp);
250: KSPSolve(ksp_A, b, x);
251: }
252: }
254: KSPDestroy(&ksp_A);
255: MatDestroy(&P);
256: MatDestroy(&A);
257: VecDestroy(&x);
258: VecDestroy(&b);
259: ISDestroy(&isu);
260: ISDestroy(&isp);
261: return 0;
262: }
264: int main(int argc, char **argv)
265: {
267: PetscInitialize(&argc, &argv, 0, help);
268: port_lsd_bfbt();
269: PetscFinalize();
270: return 0;
271: }
273: /*TEST
275: test:
276: args: -f ${DATAFILESPATH}/matrices/underworld32.gz -fc_ksp_view -fc_ksp_monitor_short -fc_ksp_type fgmres -fc_ksp_max_it 4000 -fc_ksp_diagonal_scale -fc_pc_type fieldsplit -fc_pc_fieldsplit_type SCHUR -fc_pc_fieldsplit_schur_fact_type UPPER -fc_fieldsplit_velocity_ksp_type cg -fc_fieldsplit_velocity_pc_type cholesky -fc_fieldsplit_velocity_pc_factor_mat_ordering_type nd -fc_fieldsplit_pressure_ksp_max_it 100 -fc_fieldsplit_pressure_ksp_constant_null_space -fc_fieldsplit_pressure_ksp_monitor_short -fc_fieldsplit_pressure_pc_type lsc -fc_fieldsplit_pressure_lsc_ksp_type cg -fc_fieldsplit_pressure_lsc_ksp_max_it 100 -fc_fieldsplit_pressure_lsc_ksp_constant_null_space -fc_fieldsplit_pressure_lsc_ksp_converged_reason -fc_fieldsplit_pressure_lsc_pc_type icc -test_fs {{0 1}separate output} -fc_pc_fieldsplit_off_diag_use_amat {{0 1}separate output} -fc_pc_fieldsplit_diag_use_amat {{0 1}separate output}
277: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
279: test:
280: suffix: 2
281: nsize: 4
282: args: -f ${DATAFILESPATH}/matrices/underworld32.gz -fc_ksp_view -fc_ksp_monitor_short -fc_ksp_type fgmres -fc_ksp_max_it 4000 -fc_ksp_diagonal_scale -fc_pc_type fieldsplit -fc_pc_fieldsplit_type SCHUR -fc_pc_fieldsplit_schur_fact_type UPPER -fc_fieldsplit_velocity_ksp_type cg -fc_fieldsplit_velocity_ksp_rtol 1.0e-6 -fc_fieldsplit_velocity_pc_type bjacobi -fc_fieldsplit_velocity_sub_pc_type cholesky -fc_fieldsplit_velocity_sub_pc_factor_mat_ordering_type nd -fc_fieldsplit_pressure_ksp_type fgmres -fc_fieldsplit_pressure_ksp_constant_null_space -fc_fieldsplit_pressure_ksp_monitor_short -fc_fieldsplit_pressure_pc_type lsc -fc_fieldsplit_pressure_lsc_ksp_type cg -fc_fieldsplit_pressure_lsc_ksp_rtol 1.0e-2 -fc_fieldsplit_pressure_lsc_ksp_constant_null_space -fc_fieldsplit_pressure_lsc_ksp_converged_reason -fc_fieldsplit_pressure_lsc_pc_type bjacobi -fc_fieldsplit_pressure_lsc_sub_pc_type icc -test_fs {{0 1}separate output} -fc_pc_fieldsplit_off_diag_use_amat {{0 1}separate output} -fc_pc_fieldsplit_diag_use_amat {{0 1}separate output}
283: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
285: test:
286: suffix: 3
287: nsize: 2
288: args: -f ${DATAFILESPATH}/matrices/underworld32.gz -fc_ksp_view_pre -fc_pc_type lu
289: requires: datafilespath mumps double !complex !defined(PETSC_USE_64BIT_INDICES)
291: TEST*/