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;

 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: }

 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;

 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: }

 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;

 83:   PetscFunctionBeginUser;
 84:   /* fetch test matrices and vectors */
 85:   PetscCall(LSCLoadTestOperators(&Auu, &Aup, &Apu, &App, &f, &h));

 87:   /* build the mat-nest */
 88:   PetscCall(VecGetSize(f, &nu));
 89:   PetscCall(VecGetSize(h, &np));

 91:   PetscCall(VecGetLocalSize(f, &lnu));
 92:   PetscCall(VecGetLocalSize(h, &lnp));

 94:   PetscCall(VecGetOwnershipRange(f, &start_u, &end_u));
 95:   PetscCall(VecGetOwnershipRange(h, &start_p, &end_p));

 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));

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));

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));

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:   }

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));

138:   *_isu = is_u;
139:   *_isp = is_p;
140:   *_A   = A;
141:   *_x   = x;
142:   *_b   = b;
143:   PetscFunctionReturn(PETSC_SUCCESS);
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;

155:   PetscFunctionBeginUser;
156:   PetscCall(LoadTestMatrices(&A, &x, &b, &isu, &isp));
157:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_fs", &test_fs, NULL));
158:   if (!test_fs) {
159:     PetscCall(PetscObjectReference((PetscObject)A));
160:     P = A;
161:   } else {
162:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &P));
163:   }
164:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_A));
165:   PetscCall(KSPSetOptionsPrefix(ksp_A, "fc_"));
166:   PetscCall(KSPSetOperators(ksp_A, A, P));

168:   PetscCall(KSPSetFromOptions(ksp_A));
169:   PetscCall(KSPGetPC(ksp_A, &pc_A));
170:   PetscCall(PetscObjectTypeCompare((PetscObject)pc_A, PCLU, &test_fs));
171:   if (test_fs) {
172:     PetscCall(MatDestroy(&P));
173:     PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &P));
174:     PetscCall(KSPSetOperators(ksp_A, A, P));
175:     PetscCall(KSPSetFromOptions(ksp_A));
176:     PetscCall(KSPSolve(ksp_A, b, x));
177:   } else {
178:     PetscCall(PCFieldSplitSetBlockSize(pc_A, 2));
179:     PetscCall(PCFieldSplitSetIS(pc_A, "velocity", isu));
180:     PetscCall(PCFieldSplitSetIS(pc_A, "pressure", isp));
181:     PetscCall(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:       PetscCall(MatCreateSubMatrix(A, isu, isu, MAT_INITIAL_MATRIX, &A11));
194:       PetscCall(MatCreateSubMatrix(A, isp, isp, MAT_INITIAL_MATRIX, &A22));

196:       PetscCall(MatCreateVecs(A11, &uvec, NULL));
197:       PetscCall(MatCreateVecs(A22, &pvec, NULL));

199:       /* perform the scatter from x -> (u,p) */
200:       PetscCall(VecScatterCreate(x, isu, uvec, NULL, &uscat));
201:       PetscCall(VecScatterBegin(uscat, x, uvec, INSERT_VALUES, SCATTER_FORWARD));
202:       PetscCall(VecScatterEnd(uscat, x, uvec, INSERT_VALUES, SCATTER_FORWARD));

204:       PetscCall(VecScatterCreate(x, isp, pvec, NULL, &pscat));
205:       PetscCall(VecScatterBegin(pscat, x, pvec, INSERT_VALUES, SCATTER_FORWARD));
206:       PetscCall(VecScatterEnd(pscat, x, pvec, INSERT_VALUES, SCATTER_FORWARD));

208:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "-- vector vector values --\n"));
209:       PetscCall(VecMin(uvec, &loc, &max));
210:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Min(u)  = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
211:       PetscCall(VecMax(uvec, &loc, &max));
212:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Max(u)  = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
213:       PetscCall(VecNorm(uvec, NORM_2, &norm));
214:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Norm(u) = %1.6f \n", (double)norm));
215:       PetscCall(VecSum(uvec, &sum));
216:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Sum(u)  = %1.6f \n", (double)PetscRealPart(sum)));

218:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "-- pressure vector values --\n"));
219:       PetscCall(VecMin(pvec, &loc, &max));
220:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Min(p)  = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
221:       PetscCall(VecMax(pvec, &loc, &max));
222:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Max(p)  = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
223:       PetscCall(VecNorm(pvec, NORM_2, &norm));
224:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Norm(p) = %1.6f \n", (double)norm));
225:       PetscCall(VecSum(pvec, &sum));
226:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Sum(p)  = %1.6f \n", (double)PetscRealPart(sum)));

228:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "-- Full vector values --\n"));
229:       PetscCall(VecMin(x, &loc, &max));
230:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Min(u,p)  = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
231:       PetscCall(VecMax(x, &loc, &max));
232:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Max(u,p)  = %1.6f [loc=%" PetscInt_FMT "]\n", (double)max, loc));
233:       PetscCall(VecNorm(x, NORM_2, &norm));
234:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Norm(u,p) = %1.6f \n", (double)norm));
235:       PetscCall(VecSum(x, &sum));
236:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  Sum(u,p)  = %1.6f \n", (double)PetscRealPart(sum)));

238:       PetscCall(VecScatterDestroy(&uscat));
239:       PetscCall(VecScatterDestroy(&pscat));
240:       PetscCall(VecDestroy(&uvec));
241:       PetscCall(VecDestroy(&pvec));
242:       PetscCall(MatDestroy(&A11));
243:       PetscCall(MatDestroy(&A22));
244:     }

246:     /* test second solve by changing the mat associated to the MATNEST blocks */
247:     {
248:       PetscCall(replace_submats(A, isu, isp));
249:       PetscCall(replace_submats(P, isu, isp));
250:       PetscCall(KSPSolve(ksp_A, b, x));
251:     }
252:   }

254:   PetscCall(KSPDestroy(&ksp_A));
255:   PetscCall(MatDestroy(&P));
256:   PetscCall(MatDestroy(&A));
257:   PetscCall(VecDestroy(&x));
258:   PetscCall(VecDestroy(&b));
259:   PetscCall(ISDestroy(&isu));
260:   PetscCall(ISDestroy(&isp));
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }

264: int main(int argc, char **argv)
265: {
266:   PetscFunctionBeginUser;
267:   PetscCall(PetscInitialize(&argc, &argv, 0, help));
268:   PetscCall(port_lsd_bfbt());
269:   PetscCall(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*/