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*/