1 static const char help[] = "Solves a Q1-P0 Stokes problem from Underworld.\n\
2 \n\
3 You can obtain a sample matrix from https://web.cels.anl.gov/projects/petsc/download/Datafiles/matrices/underworld32.gz\n\
4 and run with -f underworld32.gz\n\n";
5
6 #include <petscksp.h>
7 #include <petscdmda.h>
8
replace_submats(Mat A,IS isu,IS isp)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;
14
15 PetscFunctionBeginUser;
16 PetscCall(MatCreateSubMatrix(A, isu, isu, MAT_INITIAL_MATRIX, &A11));
17 PetscCall(MatCreateSubMatrix(A, isu, isp, MAT_INITIAL_MATRIX, &A12));
18 PetscCall(MatCreateSubMatrix(A, isp, isu, MAT_INITIAL_MATRIX, &A21));
19 PetscCall(MatCreateSubMatrix(A, isp, isp, MAT_INITIAL_MATRIX, &A22));
20 PetscCall(MatDuplicate(A11, MAT_COPY_VALUES, &nA11));
21 PetscCall(MatDuplicate(A12, MAT_COPY_VALUES, &nA12));
22 PetscCall(MatDuplicate(A21, MAT_COPY_VALUES, &nA21));
23 PetscCall(MatDuplicate(A22, MAT_COPY_VALUES, &nA22));
24 PetscCall(MatGetOptionsPrefix(A11, &prefix));
25 PetscCall(MatSetOptionsPrefix(nA11, prefix));
26 PetscCall(MatGetOptionsPrefix(A22, &prefix));
27 PetscCall(MatSetOptionsPrefix(nA22, prefix));
28 PetscCall(MatNestSetSubMat(A, 0, 0, nA11));
29 PetscCall(MatNestSetSubMat(A, 0, 1, nA12));
30 PetscCall(MatNestSetSubMat(A, 1, 0, nA21));
31 PetscCall(MatNestSetSubMat(A, 1, 1, nA22));
32 PetscCall(MatDestroy(&A11));
33 PetscCall(MatDestroy(&A12));
34 PetscCall(MatDestroy(&A21));
35 PetscCall(MatDestroy(&A22));
36 PetscCall(MatDestroy(&nA11));
37 PetscCall(MatDestroy(&nA12));
38 PetscCall(MatDestroy(&nA21));
39 PetscCall(MatDestroy(&nA22));
40 PetscFunctionReturn(PETSC_SUCCESS);
41 }
42
LSCLoadTestOperators(Mat * A11,Mat * A12,Mat * A21,Mat * A22,Vec * b1,Vec * b2)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;
48
49 PetscFunctionBeginUser;
50 PetscCall(MatCreate(PETSC_COMM_WORLD, A11));
51 PetscCall(MatCreate(PETSC_COMM_WORLD, A12));
52 PetscCall(MatCreate(PETSC_COMM_WORLD, A21));
53 PetscCall(MatCreate(PETSC_COMM_WORLD, A22));
54 PetscCall(MatSetOptionsPrefix(*A11, "a11_"));
55 PetscCall(MatSetOptionsPrefix(*A22, "a22_"));
56 PetscCall(MatSetFromOptions(*A11));
57 PetscCall(MatSetFromOptions(*A22));
58 PetscCall(VecCreate(PETSC_COMM_WORLD, b1));
59 PetscCall(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 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
62 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must provide a matrix file with -f");
63 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer));
64 PetscCall(MatLoad(*A11, viewer));
65 PetscCall(MatLoad(*A12, viewer));
66 PetscCall(MatLoad(*A21, viewer));
67 PetscCall(MatLoad(*A22, viewer));
68 PetscCall(VecLoad(*b1, viewer));
69 PetscCall(VecLoad(*b2, viewer));
70 PetscCall(PetscViewerDestroy(&viewer));
71 PetscFunctionReturn(PETSC_SUCCESS);
72 }
73
LoadTestMatrices(Mat * _A,Vec * _x,Vec * _b,IS * _isu,IS * _isp)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;
82
83 PetscFunctionBeginUser;
84 /* fetch test matrices and vectors */
85 PetscCall(LSCLoadTestOperators(&Auu, &Aup, &Apu, &App, &f, &h));
86
87 /* build the mat-nest */
88 PetscCall(VecGetSize(f, &nu));
89 PetscCall(VecGetSize(h, &np));
90
91 PetscCall(VecGetLocalSize(f, &lnu));
92 PetscCall(VecGetLocalSize(h, &lnp));
93
94 PetscCall(VecGetOwnershipRange(f, &start_u, &end_u));
95 PetscCall(VecGetOwnershipRange(h, &start_p, &end_p));
96
97 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
98 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] lnu = %" PetscInt_FMT " | lnp = %" PetscInt_FMT " \n", rank, lnu, lnp));
99 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] s_u = %" PetscInt_FMT " | e_u = %" PetscInt_FMT " \n", rank, start_u, end_u));
100 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] s_p = %" PetscInt_FMT " | e_p = %" PetscInt_FMT " \n", rank, start_p, end_p));
101 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] is_u (offset) = %" PetscInt_FMT " \n", rank, start_u + start_p));
102 PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] is_p (offset) = %" PetscInt_FMT " \n", rank, start_u + start_p + lnu));
103 PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
104
105 PetscCall(ISCreateStride(PETSC_COMM_WORLD, lnu, start_u + start_p, 1, &is_u));
106 PetscCall(ISCreateStride(PETSC_COMM_WORLD, lnp, start_u + start_p + lnu, 1, &is_p));
107
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 PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, bis, 2, bis, &bA[0][0], &A));
115 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
116 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
117
118 /* Pull f,h into b */
119 PetscCall(MatCreateVecs(A, &b, &x));
120 bX[0] = f;
121 bX[1] = h;
122 PetscCall(PetscMalloc1(2, &vscat));
123 for (i = 0; i < 2; i++) {
124 PetscCall(VecScatterCreate(b, bis[i], bX[i], NULL, &vscat[i]));
125 PetscCall(VecScatterBegin(vscat[i], bX[i], b, INSERT_VALUES, SCATTER_REVERSE));
126 PetscCall(VecScatterEnd(vscat[i], bX[i], b, INSERT_VALUES, SCATTER_REVERSE));
127 PetscCall(VecScatterDestroy(&vscat[i]));
128 }
129
130 PetscCall(PetscFree(vscat));
131 PetscCall(MatDestroy(&Auu));
132 PetscCall(MatDestroy(&Aup));
133 PetscCall(MatDestroy(&Apu));
134 PetscCall(MatDestroy(&App));
135 PetscCall(VecDestroy(&f));
136 PetscCall(VecDestroy(&h));
137
138 *_isu = is_u;
139 *_isp = is_p;
140 *_A = A;
141 *_x = x;
142 *_b = b;
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
port_lsd_bfbt(void)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, set_symmetry = PETSC_FALSE, test_sbaij = PETSC_FALSE;
154
155 PetscFunctionBeginUser;
156 PetscCall(LoadTestMatrices(&A, &x, &b, &isu, &isp));
157 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_fs", &test_fs, NULL));
158 PetscCall(PetscOptionsGetBool(NULL, NULL, "-set_symmetry", &set_symmetry, NULL));
159 PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_sbaij", &test_sbaij, NULL));
160 if (!test_fs) {
161 PetscCall(PetscObjectReference((PetscObject)A));
162 P = A;
163 } else {
164 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &P));
165 }
166 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_A));
167 PetscCall(KSPSetOptionsPrefix(ksp_A, "fc_"));
168 PetscCall(KSPSetOperators(ksp_A, A, P));
169
170 PetscCall(KSPSetFromOptions(ksp_A));
171 PetscCall(KSPGetPC(ksp_A, &pc_A));
172 PetscCall(PetscObjectTypeCompare((PetscObject)pc_A, PCFIELDSPLIT, &test_fs));
173 if (!test_fs) {
174 PetscCall(MatDestroy(&P));
175 PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &P));
176 if (set_symmetry || test_sbaij) {
177 PetscCall(MatSetOption(P, MAT_SYMMETRIC, PETSC_TRUE));
178 if (test_sbaij) PetscCall(MatConvert(P, MATSBAIJ, MAT_INPLACE_MATRIX, &P));
179 }
180 PetscCall(KSPSetOperators(ksp_A, A, P));
181 PetscCall(KSPSetFromOptions(ksp_A));
182 PetscCall(KSPSolve(ksp_A, b, x));
183 } else {
184 PetscCall(PCFieldSplitSetBlockSize(pc_A, 2));
185 PetscCall(PCFieldSplitSetIS(pc_A, "velocity", isu));
186 PetscCall(PCFieldSplitSetIS(pc_A, "pressure", isp));
187 PetscCall(KSPSolve(ksp_A, b, x));
188
189 /* Pull u,p out of x */
190 {
191 PetscInt loc;
192 PetscReal max, norm;
193 PetscScalar sum;
194 Vec uvec, pvec;
195 VecScatter uscat, pscat;
196 Mat A11, A22;
197
198 /* grab matrices and create the compatible u,p vectors */
199 PetscCall(MatCreateSubMatrix(A, isu, isu, MAT_INITIAL_MATRIX, &A11));
200 PetscCall(MatCreateSubMatrix(A, isp, isp, MAT_INITIAL_MATRIX, &A22));
201
202 PetscCall(MatCreateVecs(A11, &uvec, NULL));
203 PetscCall(MatCreateVecs(A22, &pvec, NULL));
204
205 /* perform the scatter from x -> (u,p) */
206 PetscCall(VecScatterCreate(x, isu, uvec, NULL, &uscat));
207 PetscCall(VecScatterBegin(uscat, x, uvec, INSERT_VALUES, SCATTER_FORWARD));
208 PetscCall(VecScatterEnd(uscat, x, uvec, INSERT_VALUES, SCATTER_FORWARD));
209
210 PetscCall(VecScatterCreate(x, isp, pvec, NULL, &pscat));
211 PetscCall(VecScatterBegin(pscat, x, pvec, INSERT_VALUES, SCATTER_FORWARD));
212 PetscCall(VecScatterEnd(pscat, x, pvec, INSERT_VALUES, SCATTER_FORWARD));
213
214 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "-- vector vector values --\n"));
215 PetscCall(VecMin(uvec, &loc, &max));
216 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Min(u) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
217 PetscCall(VecMax(uvec, &loc, &max));
218 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Max(u) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
219 PetscCall(VecNorm(uvec, NORM_2, &norm));
220 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Norm(u) = %1.6f \n", (double)norm));
221 PetscCall(VecSum(uvec, &sum));
222 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Sum(u) = %1.6f \n", (double)PetscRealPart(sum)));
223
224 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "-- pressure vector values --\n"));
225 PetscCall(VecMin(pvec, &loc, &max));
226 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Min(p) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
227 PetscCall(VecMax(pvec, &loc, &max));
228 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Max(p) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
229 PetscCall(VecNorm(pvec, NORM_2, &norm));
230 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Norm(p) = %1.6f \n", (double)norm));
231 PetscCall(VecSum(pvec, &sum));
232 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Sum(p) = %1.6f \n", (double)PetscRealPart(sum)));
233
234 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "-- Full vector values --\n"));
235 PetscCall(VecMin(x, &loc, &max));
236 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Min(u,p) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
237 PetscCall(VecMax(x, &loc, &max));
238 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Max(u,p) = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
239 PetscCall(VecNorm(x, NORM_2, &norm));
240 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Norm(u,p) = %1.6f \n", (double)norm));
241 PetscCall(VecSum(x, &sum));
242 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Sum(u,p) = %1.6f \n", (double)PetscRealPart(sum)));
243
244 PetscCall(VecScatterDestroy(&uscat));
245 PetscCall(VecScatterDestroy(&pscat));
246 PetscCall(VecDestroy(&uvec));
247 PetscCall(VecDestroy(&pvec));
248 PetscCall(MatDestroy(&A11));
249 PetscCall(MatDestroy(&A22));
250 }
251
252 /* test second solve by changing the mat associated to the MATNEST blocks */
253 {
254 PetscCall(replace_submats(A, isu, isp));
255 PetscCall(replace_submats(P, isu, isp));
256 PetscCall(KSPSolve(ksp_A, b, x));
257 }
258 }
259
260 PetscCall(KSPDestroy(&ksp_A));
261 PetscCall(MatDestroy(&P));
262 PetscCall(MatDestroy(&A));
263 PetscCall(VecDestroy(&x));
264 PetscCall(VecDestroy(&b));
265 PetscCall(ISDestroy(&isu));
266 PetscCall(ISDestroy(&isp));
267 PetscFunctionReturn(PETSC_SUCCESS);
268 }
269
main(int argc,char ** argv)270 int main(int argc, char **argv)
271 {
272 PetscFunctionBeginUser;
273 PetscCall(PetscInitialize(&argc, &argv, 0, help));
274 PetscCall(port_lsd_bfbt());
275 PetscCall(PetscFinalize());
276 return 0;
277 }
278
279 /*TEST
280
281 test:
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_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} -a11_mat_block_size 2
283 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
284
285 test:
286 suffix: 2
287 nsize: 4
288 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}
289 requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
290
291 test:
292 suffix: 3
293 nsize: 2
294 args: -f ${DATAFILESPATH}/matrices/underworld32.gz -fc_ksp_view_pre -fc_pc_type lu
295 requires: datafilespath mumps double !complex !defined(PETSC_USE_64BIT_INDICES)
296
297 testset:
298 requires: hpddm slepc datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) defined(PETSC_HAVE_DYNAMIC_LIBRARIES) defined(PETSC_USE_SHARED_LIBRARIES)
299 nsize: 4
300 args: -f ${DATAFILESPATH}/matrices/underworld32.gz -test_fs false -prefix_push fc_ -ksp_converged_reason -ksp_max_it 100 -ksp_pc_side right -pc_type hpddm -pc_hpddm_levels_1_svd_nsv 100 -pc_hpddm_levels_1_svd_threshold_relative 1e-6 -pc_hpddm_levels_1_sub_pc_type cholesky -pc_hpddm_levels_1_sub_pc_factor_shift_type inblocks -prefix_pop
301 test:
302 suffix: harmonic_overlap_1
303 filter: grep -v "WARNING! " | grep -v "There are 2 unused database options" | grep -v "Option left: name:-fc_pc_hpddm_levels_1_svd_pc_"
304 args: -prefix_push fc_ -pc_hpddm_harmonic_overlap 1 -pc_hpddm_levels_1_svd_pc_type cholesky -pc_hpddm_levels_1_svd_pc_factor_shift_type inblocks -prefix_pop -fc_pc_hpddm_levels_1_st_share_sub_ksp {{true false}shared output} -test_sbaij {{true false}shared output}
305 test:
306 suffix: harmonic_overlap_1_define_false
307 output_file: output/ex11_harmonic_overlap_1.out
308 args: -prefix_push fc_ -pc_hpddm_harmonic_overlap 1 -pc_hpddm_define_subdomains false -pc_hpddm_levels_1_svd_pc_type cholesky -pc_hpddm_levels_1_svd_pc_factor_shift_type inblocks -pc_hpddm_levels_1_svd_pc_factor_shift_type inblocks -pc_hpddm_levels_1_pc_type asm -pc_hpddm_levels_1_sub_pc_type cholesky -prefix_pop
309
310 TEST*/
311