Actual source code: ex73.c

  1: /*
  2: This example was derived from src/ksp/ksp/tutorials ex29.c

  4: Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation

  6:    -div \rho grad u = f,  0 < x,y < 1,

  8: with forcing function

 10:    f = e^{-x^2/\nu} e^{-y^2/\nu}

 12: with Dirichlet boundary conditions

 14:    u = f(x,y) for x = 0, x = 1, y = 0, y = 1

 16: or pure Neumman boundary conditions.
 17: */

 19: static char help[] = "Solves 2D inhomogeneous Laplacian. Demonstrates using PCTelescopeSetCoarseDM functionality of PCTelescope via a DMShell\n\n";

 21: #include <petscdm.h>
 22: #include <petscdmda.h>
 23: #include <petscdmshell.h>
 24: #include <petscksp.h>

 26: PetscErrorCode ComputeMatrix_ShellDA(KSP, Mat, Mat, void *);
 27: PetscErrorCode ComputeMatrix_DMDA(DM, Mat, Mat, void *);
 28: PetscErrorCode ComputeRHS_DMDA(DM, Vec, void *);
 29: PetscErrorCode DMShellCreate_ShellDA(DM, DM *);
 30: PetscErrorCode DMFieldScatter_ShellDA(DM, Vec, ScatterMode, DM, Vec);
 31: PetscErrorCode DMStateScatter_ShellDA(DM, ScatterMode, DM);

 33: typedef enum {
 34:   DIRICHLET,
 35:   NEUMANN
 36: } BCType;

 38: typedef struct {
 39:   PetscReal rho;
 40:   PetscReal nu;
 41:   BCType    bcType;
 42:   MPI_Comm  comm;
 43: } UserContext;

 45: PetscErrorCode UserContextCreate(MPI_Comm comm, UserContext **ctx)
 46: {
 47:   UserContext *user;
 48:   const char  *bcTypes[2] = {"dirichlet", "neumann"};
 49:   PetscInt     bc;

 51:   PetscFunctionBeginUser;
 52:   PetscCall(PetscCalloc1(1, &user));
 53:   user->comm = comm;
 54:   PetscOptionsBegin(comm, "", "Options for the inhomogeneous Poisson equation", "DMqq");
 55:   user->rho = 1.0;
 56:   PetscCall(PetscOptionsReal("-rho", "The conductivity", "ex29.c", user->rho, &user->rho, NULL));
 57:   user->nu = 0.1;
 58:   PetscCall(PetscOptionsReal("-nu", "The width of the Gaussian source", "ex29.c", user->nu, &user->nu, NULL));
 59:   bc = (PetscInt)DIRICHLET;
 60:   PetscCall(PetscOptionsEList("-bc_type", "Type of boundary condition", "ex29.c", bcTypes, 2, bcTypes[0], &bc, NULL));
 61:   user->bcType = (BCType)bc;
 62:   PetscOptionsEnd();
 63:   *ctx = user;
 64:   PetscFunctionReturn(PETSC_SUCCESS);
 65: }

 67: PetscErrorCode CommCoarsen(MPI_Comm comm, PetscInt number, PetscSubcomm *p)
 68: {
 69:   PetscSubcomm psubcomm;
 70:   PetscFunctionBeginUser;
 71:   PetscCall(PetscSubcommCreate(comm, &psubcomm));
 72:   PetscCall(PetscSubcommSetNumber(psubcomm, number));
 73:   PetscCall(PetscSubcommSetType(psubcomm, PETSC_SUBCOMM_INTERLACED));
 74:   *p = psubcomm;
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }

 78: PetscErrorCode CommHierarchyCreate(MPI_Comm comm, PetscInt n, PetscInt number[], PetscSubcomm pscommlist[])
 79: {
 80:   PetscInt  k;
 81:   PetscBool view_hierarchy = PETSC_FALSE;

 83:   PetscFunctionBeginUser;
 84:   for (k = 0; k < n; k++) pscommlist[k] = NULL;

 86:   if (n < 1) PetscFunctionReturn(PETSC_SUCCESS);

 88:   PetscCall(CommCoarsen(comm, number[n - 1], &pscommlist[n - 1]));
 89:   for (k = n - 2; k >= 0; k--) {
 90:     MPI_Comm comm_k = PetscSubcommChild(pscommlist[k + 1]);
 91:     if (pscommlist[k + 1]->color == 0) PetscCall(CommCoarsen(comm_k, number[k], &pscommlist[k]));
 92:   }

 94:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_hierarchy", &view_hierarchy, NULL));
 95:   if (view_hierarchy) {
 96:     PetscMPIInt size;

 98:     PetscCallMPI(MPI_Comm_size(comm, &size));
 99:     PetscCall(PetscPrintf(comm, "level[%" PetscInt_FMT "] size %d\n", n, (int)size));
100:     for (k = n - 1; k >= 0; k--) {
101:       if (pscommlist[k]) {
102:         MPI_Comm comm_k = PetscSubcommChild(pscommlist[k]);

104:         if (pscommlist[k]->color == 0) {
105:           PetscCallMPI(MPI_Comm_size(comm_k, &size));
106:           PetscCall(PetscPrintf(comm_k, "level[%" PetscInt_FMT "] size %d\n", k, (int)size));
107:         }
108:       }
109:     }
110:   }
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }

114: /* taken from src/ksp/pc/impls/telescope/telescope_dmda.c */
115: static PetscErrorCode _DMDADetermineRankFromGlobalIJ_2d(PetscInt i, PetscInt j, PetscInt Mp, PetscInt Np, PetscInt start_i[], PetscInt start_j[], PetscInt span_i[], PetscInt span_j[], PetscMPIInt *_pi, PetscMPIInt *_pj, PetscMPIInt *rank_re)
116: {
117:   PetscInt pi, pj, n;

119:   PetscFunctionBeginUser;
120:   *rank_re = -1;
121:   pi = pj = -1;
122:   if (_pi) {
123:     for (n = 0; n < Mp; n++) {
124:       if ((i >= start_i[n]) && (i < start_i[n] + span_i[n])) {
125:         pi = n;
126:         break;
127:       }
128:     }
129:     PetscCheck(pi != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmda-ij] pi cannot be determined : range %" PetscInt_FMT ", val %" PetscInt_FMT, Mp, i);
130:     *_pi = (PetscMPIInt)pi;
131:   }

133:   if (_pj) {
134:     for (n = 0; n < Np; n++) {
135:       if ((j >= start_j[n]) && (j < start_j[n] + span_j[n])) {
136:         pj = n;
137:         break;
138:       }
139:     }
140:     PetscCheck(pj != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmda-ij] pj cannot be determined : range %" PetscInt_FMT ", val %" PetscInt_FMT, Np, j);
141:     *_pj = (PetscMPIInt)pj;
142:   }

144:   *rank_re = (PetscMPIInt)(pi + pj * Mp);
145:   PetscFunctionReturn(PETSC_SUCCESS);
146: }

148: /* taken from src/ksp/pc/impls/telescope/telescope_dmda.c */
149: static PetscErrorCode _DMDADetermineGlobalS0_2d(PetscMPIInt rank_re, PetscInt Mp_re, PetscInt Np_re, PetscInt range_i_re[], PetscInt range_j_re[], PetscInt *s0)
150: {
151:   PetscInt    i, j, start_IJ = 0;
152:   PetscMPIInt rank_ij;

154:   PetscFunctionBeginUser;
155:   *s0 = -1;
156:   for (j = 0; j < Np_re; j++) {
157:     for (i = 0; i < Mp_re; i++) {
158:       rank_ij = (PetscMPIInt)(i + j * Mp_re);
159:       if (rank_ij < rank_re) start_IJ += range_i_re[i] * range_j_re[j];
160:     }
161:   }
162:   *s0 = start_IJ;
163:   PetscFunctionReturn(PETSC_SUCCESS);
164: }

166: /* adapted from src/ksp/pc/impls/telescope/telescope_dmda.c */
167: static PetscErrorCode DMDACreatePermutation_2d(DM dmrepart, DM dmf, Mat *mat)
168: {
169:   PetscInt        k, sum, Mp_re = 0, Np_re = 0;
170:   PetscInt        nx, ny, sr, er, Mr, ndof;
171:   PetscInt        i, j, location, startI[2], endI[2], lenI[2];
172:   const PetscInt *_range_i_re = NULL, *_range_j_re = NULL;
173:   PetscInt       *range_i_re, *range_j_re;
174:   PetscInt       *start_i_re, *start_j_re;
175:   MPI_Comm        comm;
176:   Vec             V;
177:   Mat             Pscalar;

179:   PetscFunctionBeginUser;
180:   PetscCall(PetscInfo(dmf, "setting up the permutation matrix (DMDA-2D)\n"));
181:   PetscCall(PetscObjectGetComm((PetscObject)dmf, &comm));

183:   _range_i_re = _range_j_re = NULL;
184:   /* Create DMDA on the child communicator */
185:   if (dmrepart) {
186:     PetscCall(DMDAGetInfo(dmrepart, NULL, NULL, NULL, NULL, &Mp_re, &Np_re, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
187:     PetscCall(DMDAGetOwnershipRanges(dmrepart, &_range_i_re, &_range_j_re, NULL));
188:   }

190:   /* note - assume rank 0 always participates */
191:   PetscCallMPI(MPI_Bcast(&Mp_re, 1, MPIU_INT, 0, comm));
192:   PetscCallMPI(MPI_Bcast(&Np_re, 1, MPIU_INT, 0, comm));

194:   PetscCall(PetscCalloc1(Mp_re, &range_i_re));
195:   PetscCall(PetscCalloc1(Np_re, &range_j_re));

197:   if (_range_i_re) PetscCall(PetscArraycpy(range_i_re, _range_i_re, Mp_re));
198:   if (_range_j_re) PetscCall(PetscArraycpy(range_j_re, _range_j_re, Np_re));

200:   PetscCallMPI(MPI_Bcast(range_i_re, Mp_re, MPIU_INT, 0, comm));
201:   PetscCallMPI(MPI_Bcast(range_j_re, Np_re, MPIU_INT, 0, comm));

203:   PetscCall(PetscMalloc1(Mp_re, &start_i_re));
204:   PetscCall(PetscMalloc1(Np_re, &start_j_re));

206:   sum = 0;
207:   for (k = 0; k < Mp_re; k++) {
208:     start_i_re[k] = sum;
209:     sum += range_i_re[k];
210:   }

212:   sum = 0;
213:   for (k = 0; k < Np_re; k++) {
214:     start_j_re[k] = sum;
215:     sum += range_j_re[k];
216:   }

218:   /* Create permutation */
219:   PetscCall(DMDAGetInfo(dmf, NULL, &nx, &ny, NULL, NULL, NULL, NULL, &ndof, NULL, NULL, NULL, NULL, NULL));
220:   PetscCall(DMGetGlobalVector(dmf, &V));
221:   PetscCall(VecGetSize(V, &Mr));
222:   PetscCall(VecGetOwnershipRange(V, &sr, &er));
223:   PetscCall(DMRestoreGlobalVector(dmf, &V));
224:   sr = sr / ndof;
225:   er = er / ndof;
226:   Mr = Mr / ndof;

228:   PetscCall(MatCreate(comm, &Pscalar));
229:   PetscCall(MatSetSizes(Pscalar, (er - sr), (er - sr), Mr, Mr));
230:   PetscCall(MatSetType(Pscalar, MATAIJ));
231:   PetscCall(MatSeqAIJSetPreallocation(Pscalar, 1, NULL));
232:   PetscCall(MatMPIAIJSetPreallocation(Pscalar, 1, NULL, 1, NULL));

234:   PetscCall(DMDAGetCorners(dmf, NULL, NULL, NULL, &lenI[0], &lenI[1], NULL));
235:   PetscCall(DMDAGetCorners(dmf, &startI[0], &startI[1], NULL, &endI[0], &endI[1], NULL));
236:   endI[0] += startI[0];
237:   endI[1] += startI[1];

239:   for (j = startI[1]; j < endI[1]; j++) {
240:     for (i = startI[0]; i < endI[0]; i++) {
241:       PetscMPIInt rank_ijk_re, rank_reI[] = {0, 0};
242:       PetscInt    s0_re;
243:       PetscInt    ii, jj, local_ijk_re, mapped_ijk;
244:       PetscInt    lenI_re[] = {0, 0};

246:       location = (i - startI[0]) + (j - startI[1]) * lenI[0];
247:       PetscCall(_DMDADetermineRankFromGlobalIJ_2d(i, j, Mp_re, Np_re, start_i_re, start_j_re, range_i_re, range_j_re, &rank_reI[0], &rank_reI[1], &rank_ijk_re));

249:       PetscCall(_DMDADetermineGlobalS0_2d(rank_ijk_re, Mp_re, Np_re, range_i_re, range_j_re, &s0_re));

251:       ii = i - start_i_re[rank_reI[0]];
252:       PetscCheck(ii >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmdarepart-perm2d] index error ii");
253:       jj = j - start_j_re[rank_reI[1]];
254:       PetscCheck(jj >= 0, PETSC_COMM_SELF, PETSC_ERR_USER, "[dmdarepart-perm2d] index error jj");

256:       lenI_re[0]   = range_i_re[rank_reI[0]];
257:       lenI_re[1]   = range_j_re[rank_reI[1]];
258:       local_ijk_re = ii + jj * lenI_re[0];
259:       mapped_ijk   = s0_re + local_ijk_re;
260:       PetscCall(MatSetValue(Pscalar, sr + location, mapped_ijk, 1.0, INSERT_VALUES));
261:     }
262:   }
263:   PetscCall(MatAssemblyBegin(Pscalar, MAT_FINAL_ASSEMBLY));
264:   PetscCall(MatAssemblyEnd(Pscalar, MAT_FINAL_ASSEMBLY));

266:   *mat = Pscalar;

268:   PetscCall(PetscFree(range_i_re));
269:   PetscCall(PetscFree(range_j_re));
270:   PetscCall(PetscFree(start_i_re));
271:   PetscCall(PetscFree(start_j_re));
272:   PetscFunctionReturn(PETSC_SUCCESS);
273: }

275: /* adapted from src/ksp/pc/impls/telescope/telescope_dmda.c */
276: static PetscErrorCode PCTelescopeSetUp_dmda_scatters(DM dmf, DM dmc)
277: {
278:   Vec        xred, yred, xtmp, x, xp;
279:   VecScatter scatter;
280:   IS         isin;
281:   PetscInt   m, bs, st, ed;
282:   MPI_Comm   comm;
283:   VecType    vectype;

285:   PetscFunctionBeginUser;
286:   PetscCall(PetscObjectGetComm((PetscObject)dmf, &comm));
287:   PetscCall(DMCreateGlobalVector(dmf, &x));
288:   PetscCall(VecGetBlockSize(x, &bs));
289:   PetscCall(VecGetType(x, &vectype));

291:   /* cannot use VecDuplicate as xp is already composed with dmf */
292:   /*PetscCall(VecDuplicate(x,&xp));*/
293:   {
294:     PetscInt m, M;

296:     PetscCall(VecGetSize(x, &M));
297:     PetscCall(VecGetLocalSize(x, &m));
298:     PetscCall(VecCreate(comm, &xp));
299:     PetscCall(VecSetSizes(xp, m, M));
300:     PetscCall(VecSetBlockSize(xp, bs));
301:     PetscCall(VecSetType(xp, vectype));
302:   }

304:   m    = 0;
305:   xred = NULL;
306:   yred = NULL;
307:   if (dmc) {
308:     PetscCall(DMCreateGlobalVector(dmc, &xred));
309:     PetscCall(VecDuplicate(xred, &yred));
310:     PetscCall(VecGetOwnershipRange(xred, &st, &ed));
311:     PetscCall(ISCreateStride(comm, ed - st, st, 1, &isin));
312:     PetscCall(VecGetLocalSize(xred, &m));
313:   } else {
314:     PetscCall(VecGetOwnershipRange(x, &st, &ed));
315:     PetscCall(ISCreateStride(comm, 0, st, 1, &isin));
316:   }
317:   PetscCall(ISSetBlockSize(isin, bs));
318:   PetscCall(VecCreate(comm, &xtmp));
319:   PetscCall(VecSetSizes(xtmp, m, PETSC_DECIDE));
320:   PetscCall(VecSetBlockSize(xtmp, bs));
321:   PetscCall(VecSetType(xtmp, vectype));
322:   PetscCall(VecScatterCreate(x, isin, xtmp, NULL, &scatter));

324:   PetscCall(PetscObjectCompose((PetscObject)dmf, "isin", (PetscObject)isin));
325:   PetscCall(PetscObjectCompose((PetscObject)dmf, "scatter", (PetscObject)scatter));
326:   PetscCall(PetscObjectCompose((PetscObject)dmf, "xtmp", (PetscObject)xtmp));
327:   PetscCall(PetscObjectCompose((PetscObject)dmf, "xp", (PetscObject)xp));

329:   PetscCall(VecDestroy(&xred));
330:   PetscCall(VecDestroy(&yred));
331:   PetscCall(VecDestroy(&x));
332:   PetscFunctionReturn(PETSC_SUCCESS);
333: }

335: PetscErrorCode DMCreateMatrix_ShellDA(DM dm, Mat *A)
336: {
337:   DM           da;
338:   MPI_Comm     comm;
339:   PetscMPIInt  size;
340:   UserContext *ctx = NULL;
341:   PetscInt     M, N;

343:   PetscFunctionBeginUser;
344:   PetscCall(DMShellGetContext(dm, &da));
345:   PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
346:   PetscCallMPI(MPI_Comm_size(comm, &size));
347:   PetscCall(DMCreateMatrix(da, A));
348:   PetscCall(MatGetSize(*A, &M, &N));
349:   PetscCall(PetscPrintf(comm, "[size %" PetscInt_FMT "] DMCreateMatrix_ShellDA (%" PetscInt_FMT " x %" PetscInt_FMT ")\n", (PetscInt)size, M, N));

351:   PetscCall(DMGetApplicationContext(dm, &ctx));
352:   if (ctx->bcType == NEUMANN) {
353:     MatNullSpace nullspace = NULL;
354:     PetscCall(PetscPrintf(comm, "[size %" PetscInt_FMT "] DMCreateMatrix_ShellDA: using neumann bcs\n", (PetscInt)size));

356:     PetscCall(MatGetNullSpace(*A, &nullspace));
357:     if (!nullspace) {
358:       PetscCall(PetscPrintf(comm, "[size %" PetscInt_FMT "] DMCreateMatrix_ShellDA: operator does not have nullspace - attaching\n", (PetscInt)size));
359:       PetscCall(MatNullSpaceCreate(comm, PETSC_TRUE, 0, 0, &nullspace));
360:       PetscCall(MatSetNullSpace(*A, nullspace));
361:       PetscCall(MatNullSpaceDestroy(&nullspace));
362:     } else {
363:       PetscCall(PetscPrintf(comm, "[size %" PetscInt_FMT "] DMCreateMatrix_ShellDA: operator already has a nullspace\n", (PetscInt)size));
364:     }
365:   }
366:   PetscFunctionReturn(PETSC_SUCCESS);
367: }

369: PetscErrorCode DMCreateGlobalVector_ShellDA(DM dm, Vec *x)
370: {
371:   DM da;
372:   PetscFunctionBeginUser;
373:   PetscCall(DMShellGetContext(dm, &da));
374:   PetscCall(DMCreateGlobalVector(da, x));
375:   PetscCall(VecSetDM(*x, dm));
376:   PetscFunctionReturn(PETSC_SUCCESS);
377: }

379: PetscErrorCode DMCreateLocalVector_ShellDA(DM dm, Vec *x)
380: {
381:   DM da;
382:   PetscFunctionBeginUser;
383:   PetscCall(DMShellGetContext(dm, &da));
384:   PetscCall(DMCreateLocalVector(da, x));
385:   PetscCall(VecSetDM(*x, dm));
386:   PetscFunctionReturn(PETSC_SUCCESS);
387: }

389: PetscErrorCode DMCoarsen_ShellDA(DM dm, MPI_Comm comm, DM *dmc)
390: {
391:   PetscFunctionBeginUser;
392:   *dmc = NULL;
393:   PetscCall(DMGetCoarseDM(dm, dmc));
394:   if (!*dmc) {
395:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "The coarse DM should never be NULL. The DM hierarchy should have already been defined");
396:   } else {
397:     PetscCall(PetscObjectReference((PetscObject)(*dmc)));
398:   }
399:   PetscFunctionReturn(PETSC_SUCCESS);
400: }

402: PetscErrorCode DMCreateInterpolation_ShellDA(DM dm1, DM dm2, Mat *mat, Vec *vec)
403: {
404:   DM da1, da2;
405:   PetscFunctionBeginUser;
406:   PetscCall(DMShellGetContext(dm1, &da1));
407:   PetscCall(DMShellGetContext(dm2, &da2));
408:   PetscCall(DMCreateInterpolation(da1, da2, mat, vec));
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }

412: PetscErrorCode DMShellDASetUp_TelescopeDMScatter(DM dmf_shell, DM dmc_shell)
413: {
414:   Mat P   = NULL;
415:   DM  dmf = NULL, dmc = NULL;

417:   PetscFunctionBeginUser;
418:   PetscCall(DMShellGetContext(dmf_shell, &dmf));
419:   if (dmc_shell) PetscCall(DMShellGetContext(dmc_shell, &dmc));
420:   PetscCall(DMDACreatePermutation_2d(dmc, dmf, &P));
421:   PetscCall(PetscObjectCompose((PetscObject)dmf, "P", (PetscObject)P));
422:   PetscCall(PCTelescopeSetUp_dmda_scatters(dmf, dmc));
423:   PetscCall(PetscObjectComposeFunction((PetscObject)dmf_shell, "PCTelescopeFieldScatter", DMFieldScatter_ShellDA));
424:   PetscCall(PetscObjectComposeFunction((PetscObject)dmf_shell, "PCTelescopeStateScatter", DMStateScatter_ShellDA));
425:   PetscFunctionReturn(PETSC_SUCCESS);
426: }

428: PetscErrorCode DMShellDAFieldScatter_Forward(DM dmf, Vec x, DM dmc, Vec xc)
429: {
430:   Mat                P  = NULL;
431:   Vec                xp = NULL, xtmp = NULL;
432:   VecScatter         scatter = NULL;
433:   const PetscScalar *x_array;
434:   PetscInt           i, st, ed;

436:   PetscFunctionBeginUser;
437:   PetscCall(PetscObjectQuery((PetscObject)dmf, "P", (PetscObject *)&P));
438:   PetscCall(PetscObjectQuery((PetscObject)dmf, "xp", (PetscObject *)&xp));
439:   PetscCall(PetscObjectQuery((PetscObject)dmf, "scatter", (PetscObject *)&scatter));
440:   PetscCall(PetscObjectQuery((PetscObject)dmf, "xtmp", (PetscObject *)&xtmp));
441:   PetscCheck(P, PETSC_COMM_SELF, PETSC_ERR_USER, "Require a permutation matrix (\"P\")to be composed with the parent (fine) DM");
442:   PetscCheck(xp, PETSC_COMM_SELF, PETSC_ERR_USER, "Require \"xp\" to be composed with the parent (fine) DM");
443:   PetscCheck(scatter, PETSC_COMM_SELF, PETSC_ERR_USER, "Require \"scatter\" to be composed with the parent (fine) DM");
444:   PetscCheck(xtmp, PETSC_COMM_SELF, PETSC_ERR_USER, "Require \"xtmp\" to be composed with the parent (fine) DM");

446:   PetscCall(MatMultTranspose(P, x, xp));

448:   /* pull in vector x->xtmp */
449:   PetscCall(VecScatterBegin(scatter, xp, xtmp, INSERT_VALUES, SCATTER_FORWARD));
450:   PetscCall(VecScatterEnd(scatter, xp, xtmp, INSERT_VALUES, SCATTER_FORWARD));

452:   /* copy vector entries into xred */
453:   PetscCall(VecGetArrayRead(xtmp, &x_array));
454:   if (xc) {
455:     PetscScalar *LA_xred;
456:     PetscCall(VecGetOwnershipRange(xc, &st, &ed));

458:     PetscCall(VecGetArray(xc, &LA_xred));
459:     for (i = 0; i < ed - st; i++) LA_xred[i] = x_array[i];
460:     PetscCall(VecRestoreArray(xc, &LA_xred));
461:   }
462:   PetscCall(VecRestoreArrayRead(xtmp, &x_array));
463:   PetscFunctionReturn(PETSC_SUCCESS);
464: }

466: PetscErrorCode DMShellDAFieldScatter_Reverse(DM dmf, Vec y, DM dmc, Vec yc)
467: {
468:   Mat          P  = NULL;
469:   Vec          xp = NULL, xtmp = NULL;
470:   VecScatter   scatter = NULL;
471:   PetscScalar *array;
472:   PetscInt     i, st, ed;

474:   PetscFunctionBeginUser;
475:   PetscCall(PetscObjectQuery((PetscObject)dmf, "P", (PetscObject *)&P));
476:   PetscCall(PetscObjectQuery((PetscObject)dmf, "xp", (PetscObject *)&xp));
477:   PetscCall(PetscObjectQuery((PetscObject)dmf, "scatter", (PetscObject *)&scatter));
478:   PetscCall(PetscObjectQuery((PetscObject)dmf, "xtmp", (PetscObject *)&xtmp));

480:   PetscCheck(P, PETSC_COMM_SELF, PETSC_ERR_USER, "Require a permutation matrix (\"P\")to be composed with the parent (fine) DM");
481:   PetscCheck(xp, PETSC_COMM_SELF, PETSC_ERR_USER, "Require \"xp\" to be composed with the parent (fine) DM");
482:   PetscCheck(scatter, PETSC_COMM_SELF, PETSC_ERR_USER, "Require \"scatter\" to be composed with the parent (fine) DM");
483:   PetscCheck(xtmp, PETSC_COMM_SELF, PETSC_ERR_USER, "Require \"xtmp\" to be composed with the parent (fine) DM");

485:   /* return vector */
486:   PetscCall(VecGetArray(xtmp, &array));
487:   if (yc) {
488:     const PetscScalar *LA_yred;
489:     PetscCall(VecGetOwnershipRange(yc, &st, &ed));
490:     PetscCall(VecGetArrayRead(yc, &LA_yred));
491:     for (i = 0; i < ed - st; i++) array[i] = LA_yred[i];
492:     PetscCall(VecRestoreArrayRead(yc, &LA_yred));
493:   }
494:   PetscCall(VecRestoreArray(xtmp, &array));
495:   PetscCall(VecScatterBegin(scatter, xtmp, xp, INSERT_VALUES, SCATTER_REVERSE));
496:   PetscCall(VecScatterEnd(scatter, xtmp, xp, INSERT_VALUES, SCATTER_REVERSE));
497:   PetscCall(MatMult(P, xp, y));
498:   PetscFunctionReturn(PETSC_SUCCESS);
499: }

501: PetscErrorCode DMFieldScatter_ShellDA(DM dmf_shell, Vec x, ScatterMode mode, DM dmc_shell, Vec xc)
502: {
503:   DM dmf = NULL, dmc = NULL;

505:   PetscFunctionBeginUser;
506:   PetscCall(DMShellGetContext(dmf_shell, &dmf));
507:   if (dmc_shell) PetscCall(DMShellGetContext(dmc_shell, &dmc));
508:   if (mode == SCATTER_FORWARD) {
509:     PetscCall(DMShellDAFieldScatter_Forward(dmf, x, dmc, xc));
510:   } else if (mode == SCATTER_REVERSE) {
511:     PetscCall(DMShellDAFieldScatter_Reverse(dmf, x, dmc, xc));
512:   } else SETERRQ(PetscObjectComm((PetscObject)dmf_shell), PETSC_ERR_SUP, "Only mode = SCATTER_FORWARD, SCATTER_REVERSE supported");
513:   PetscFunctionReturn(PETSC_SUCCESS);
514: }

516: PetscErrorCode DMStateScatter_ShellDA(DM dmf_shell, ScatterMode mode, DM dmc_shell)
517: {
518:   PetscMPIInt size_f = 0, size_c = 0;
519:   PetscFunctionBeginUser;
520:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dmf_shell), &size_f));
521:   if (dmc_shell) PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dmc_shell), &size_c));
522:   if (mode == SCATTER_FORWARD) {
523:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dmf_shell), "User supplied state scatter (fine [size %d]-> coarse [size %d])\n", (int)size_f, (int)size_c));
524:   } else if (mode == SCATTER_REVERSE) {
525:   } else SETERRQ(PetscObjectComm((PetscObject)dmf_shell), PETSC_ERR_SUP, "Only mode = SCATTER_FORWARD, SCATTER_REVERSE supported");
526:   PetscFunctionReturn(PETSC_SUCCESS);
527: }

529: PetscErrorCode DMShellCreate_ShellDA(DM da, DM *dms)
530: {
531:   PetscFunctionBeginUser;
532:   if (da) {
533:     PetscCall(DMShellCreate(PetscObjectComm((PetscObject)da), dms));
534:     PetscCall(DMShellSetContext(*dms, da));
535:     PetscCall(DMShellSetCreateGlobalVector(*dms, DMCreateGlobalVector_ShellDA));
536:     PetscCall(DMShellSetCreateLocalVector(*dms, DMCreateLocalVector_ShellDA));
537:     PetscCall(DMShellSetCreateMatrix(*dms, DMCreateMatrix_ShellDA));
538:     PetscCall(DMShellSetCoarsen(*dms, DMCoarsen_ShellDA));
539:     PetscCall(DMShellSetCreateInterpolation(*dms, DMCreateInterpolation_ShellDA));
540:   } else {
541:     *dms = NULL;
542:   }
543:   PetscFunctionReturn(PETSC_SUCCESS);
544: }

546: PetscErrorCode DMDestroyShellDMDA(DM *_dm)
547: {
548:   DM dm, da = NULL;

550:   PetscFunctionBeginUser;
551:   if (!_dm) PetscFunctionReturn(PETSC_SUCCESS);
552:   dm = *_dm;
553:   if (!dm) PetscFunctionReturn(PETSC_SUCCESS);

555:   PetscCall(DMShellGetContext(dm, &da));
556:   if (da) {
557:     Vec        vec;
558:     VecScatter scatter = NULL;
559:     IS         is      = NULL;
560:     Mat        P       = NULL;

562:     PetscCall(PetscObjectQuery((PetscObject)da, "P", (PetscObject *)&P));
563:     PetscCall(MatDestroy(&P));

565:     vec = NULL;
566:     PetscCall(PetscObjectQuery((PetscObject)da, "xp", (PetscObject *)&vec));
567:     PetscCall(VecDestroy(&vec));

569:     PetscCall(PetscObjectQuery((PetscObject)da, "scatter", (PetscObject *)&scatter));
570:     PetscCall(VecScatterDestroy(&scatter));

572:     vec = NULL;
573:     PetscCall(PetscObjectQuery((PetscObject)da, "xtmp", (PetscObject *)&vec));
574:     PetscCall(VecDestroy(&vec));

576:     PetscCall(PetscObjectQuery((PetscObject)da, "isin", (PetscObject *)&is));
577:     PetscCall(ISDestroy(&is));

579:     PetscCall(DMDestroy(&da));
580:   }
581:   PetscCall(PetscObjectComposeFunction((PetscObject)dm, "PCTelescopeFieldScatter", NULL));
582:   PetscCall(PetscObjectComposeFunction((PetscObject)dm, "PCTelescopeStateScatter", NULL));
583:   PetscCall(DMDestroy(&dm));
584:   *_dm = NULL;
585:   PetscFunctionReturn(PETSC_SUCCESS);
586: }

588: PetscErrorCode HierarchyCreate_Basic(DM *dm_f, DM *dm_c, UserContext *ctx)
589: {
590:   DM          dm, dmc, dm_shell, dmc_shell;
591:   PetscMPIInt rank;

593:   PetscFunctionBeginUser;
594:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
595:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 17, 17, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &dm));
596:   PetscCall(DMSetFromOptions(dm));
597:   PetscCall(DMSetUp(dm));
598:   PetscCall(DMDASetUniformCoordinates(dm, 0, 1, 0, 1, 0, 0));
599:   PetscCall(DMDASetFieldName(dm, 0, "Pressure"));
600:   PetscCall(DMShellCreate_ShellDA(dm, &dm_shell));
601:   PetscCall(DMSetApplicationContext(dm_shell, ctx));

603:   dmc       = NULL;
604:   dmc_shell = NULL;
605:   if (rank == 0) {
606:     PetscCall(DMDACreate2d(PETSC_COMM_SELF, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 17, 17, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &dmc));
607:     PetscCall(DMSetFromOptions(dmc));
608:     PetscCall(DMSetUp(dmc));
609:     PetscCall(DMDASetUniformCoordinates(dmc, 0, 1, 0, 1, 0, 0));
610:     PetscCall(DMDASetFieldName(dmc, 0, "Pressure"));
611:     PetscCall(DMShellCreate_ShellDA(dmc, &dmc_shell));
612:     PetscCall(DMSetApplicationContext(dmc_shell, ctx));
613:   }

615:   PetscCall(DMSetCoarseDM(dm_shell, dmc_shell));
616:   PetscCall(DMShellDASetUp_TelescopeDMScatter(dm_shell, dmc_shell));

618:   *dm_f = dm_shell;
619:   *dm_c = dmc_shell;
620:   PetscFunctionReturn(PETSC_SUCCESS);
621: }

623: PetscErrorCode HierarchyCreate(PetscInt *_nd, PetscInt *_nref, MPI_Comm **_cl, DM **_dl)
624: {
625:   PetscInt      d, k, ndecomps, ncoarsen, found, nx;
626:   PetscInt      levelrefs, *number;
627:   PetscSubcomm *pscommlist;
628:   MPI_Comm     *commlist;
629:   DM           *dalist, *dmlist;
630:   PetscBool     set;

632:   PetscFunctionBeginUser;
633:   ndecomps = 1;
634:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ndecomps", &ndecomps, NULL));
635:   ncoarsen = ndecomps - 1;
636:   PetscCheck(ncoarsen >= 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "-ndecomps must be >= 1");

638:   levelrefs = 2;
639:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-level_nrefs", &levelrefs, NULL));
640:   PetscCall(PetscMalloc1(ncoarsen + 1, &number));
641:   for (k = 0; k < ncoarsen + 1; k++) number[k] = 2;
642:   found = ncoarsen;
643:   set   = PETSC_FALSE;
644:   PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-level_comm_red_factor", number, &found, &set));
645:   if (set) PetscCheck(found == ncoarsen, PETSC_COMM_WORLD, PETSC_ERR_USER, "Expected %" PetscInt_FMT " values for -level_comm_red_factor. Found %" PetscInt_FMT, ncoarsen, found);

647:   PetscCall(PetscMalloc1(ncoarsen + 1, &pscommlist));
648:   for (k = 0; k < ncoarsen + 1; k++) pscommlist[k] = NULL;

650:   PetscCall(PetscMalloc1(ndecomps, &commlist));
651:   for (k = 0; k < ndecomps; k++) commlist[k] = MPI_COMM_NULL;
652:   PetscCall(PetscMalloc1(ndecomps * levelrefs, &dalist));
653:   PetscCall(PetscMalloc1(ndecomps * levelrefs, &dmlist));
654:   for (k = 0; k < ndecomps * levelrefs; k++) {
655:     dalist[k] = NULL;
656:     dmlist[k] = NULL;
657:   }

659:   PetscCall(CommHierarchyCreate(PETSC_COMM_WORLD, ncoarsen, number, pscommlist));
660:   for (k = 0; k < ncoarsen; k++) {
661:     if (pscommlist[k]) {
662:       MPI_Comm comm_k = PetscSubcommChild(pscommlist[k]);
663:       if (pscommlist[k]->color == 0) PetscCall(PetscCommDuplicate(comm_k, &commlist[k], NULL));
664:     }
665:   }
666:   PetscCall(PetscCommDuplicate(PETSC_COMM_WORLD, &commlist[ndecomps - 1], NULL));

668:   for (k = 0; k < ncoarsen; k++) {
669:     if (pscommlist[k]) PetscCall(PetscSubcommDestroy(&pscommlist[k]));
670:   }

672:   nx = 17;
673:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &nx, NULL));
674:   for (d = 0; d < ndecomps; d++) {
675:     DM   dmroot = NULL;
676:     char name[PETSC_MAX_PATH_LEN];

678:     if (commlist[d] != MPI_COMM_NULL) {
679:       PetscCall(DMDACreate2d(commlist[d], DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, nx, nx, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &dmroot));
680:       PetscCall(DMSetUp(dmroot));
681:       PetscCall(DMDASetUniformCoordinates(dmroot, 0, 1, 0, 1, 0, 0));
682:       PetscCall(DMDASetFieldName(dmroot, 0, "Pressure"));
683:       PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "root-decomp-%" PetscInt_FMT, d));
684:       PetscCall(PetscObjectSetName((PetscObject)dmroot, name));
685:       /*PetscCall(DMView(dmroot,PETSC_VIEWER_STDOUT_(commlist[d])));*/
686:     }

688:     dalist[d * levelrefs + 0] = dmroot;
689:     for (k = 1; k < levelrefs; k++) {
690:       DM dmref = NULL;

692:       if (commlist[d] != MPI_COMM_NULL) {
693:         PetscCall(DMRefine(dalist[d * levelrefs + (k - 1)], MPI_COMM_NULL, &dmref));
694:         PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "ref%" PetscInt_FMT "-decomp-%" PetscInt_FMT, k, d));
695:         PetscCall(PetscObjectSetName((PetscObject)dmref, name));
696:         PetscCall(DMDAGetInfo(dmref, NULL, &nx, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
697:         /*PetscCall(DMView(dmref,PETSC_VIEWER_STDOUT_(commlist[d])));*/
698:       }
699:       dalist[d * levelrefs + k] = dmref;
700:     }
701:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, &nx, 1, MPIU_INT, MPI_MAX, PETSC_COMM_WORLD));
702:   }

704:   /* create the hierarchy of DMShell's */
705:   for (d = 0; d < ndecomps; d++) {
706:     char name[PETSC_MAX_PATH_LEN];

708:     UserContext *ctx = NULL;
709:     if (commlist[d] != MPI_COMM_NULL) {
710:       PetscCall(UserContextCreate(commlist[d], &ctx));
711:       for (k = 0; k < levelrefs; k++) {
712:         PetscCall(DMShellCreate_ShellDA(dalist[d * levelrefs + k], &dmlist[d * levelrefs + k]));
713:         PetscCall(DMSetApplicationContext(dmlist[d * levelrefs + k], ctx));
714:         PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "level%" PetscInt_FMT "-decomp-%" PetscInt_FMT, k, d));
715:         PetscCall(PetscObjectSetName((PetscObject)dmlist[d * levelrefs + k], name));
716:       }
717:     }
718:   }

720:   /* set all the coarse DMs */
721:   for (k = 1; k < ndecomps * levelrefs; k++) { /* skip first DM as it doesn't have a coarse representation */
722:     DM dmfine = NULL, dmcoarse = NULL;

724:     dmfine   = dmlist[k];
725:     dmcoarse = dmlist[k - 1];
726:     if (dmfine) PetscCall(DMSetCoarseDM(dmfine, dmcoarse));
727:   }

729:   /* do special setup on the fine DM coupling different decompositions */
730:   for (d = 1; d < ndecomps; d++) { /* skip first decomposition as it doesn't have a coarse representation */
731:     DM dmfine = NULL, dmcoarse = NULL;

733:     dmfine   = dmlist[d * levelrefs + 0];
734:     dmcoarse = dmlist[(d - 1) * levelrefs + (levelrefs - 1)];
735:     if (dmfine) PetscCall(DMShellDASetUp_TelescopeDMScatter(dmfine, dmcoarse));
736:   }

738:   PetscCall(PetscFree(number));
739:   for (k = 0; k < ncoarsen; k++) PetscCall(PetscSubcommDestroy(&pscommlist[k]));
740:   PetscCall(PetscFree(pscommlist));

742:   if (_nd) *_nd = ndecomps;
743:   if (_nref) *_nref = levelrefs;
744:   if (_cl) {
745:     *_cl = commlist;
746:   } else {
747:     for (k = 0; k < ndecomps; k++) {
748:       if (commlist[k] != MPI_COMM_NULL) PetscCall(PetscCommDestroy(&commlist[k]));
749:     }
750:     PetscCall(PetscFree(commlist));
751:   }
752:   if (_dl) {
753:     *_dl = dmlist;
754:     PetscCall(PetscFree(dalist));
755:   } else {
756:     for (k = 0; k < ndecomps * levelrefs; k++) {
757:       PetscCall(DMDestroy(&dmlist[k]));
758:       PetscCall(DMDestroy(&dalist[k]));
759:     }
760:     PetscCall(PetscFree(dmlist));
761:     PetscCall(PetscFree(dalist));
762:   }
763:   PetscFunctionReturn(PETSC_SUCCESS);
764: }

766: PetscErrorCode test_hierarchy(void)
767: {
768:   PetscInt  d, k, nd, nref;
769:   MPI_Comm *comms;
770:   DM       *dms;

772:   PetscFunctionBeginUser;
773:   PetscCall(HierarchyCreate(&nd, &nref, &comms, &dms));

775:   /* destroy user context */
776:   for (d = 0; d < nd; d++) {
777:     DM first = dms[d * nref + 0];

779:     if (first) {
780:       UserContext *ctx = NULL;

782:       PetscCall(DMGetApplicationContext(first, &ctx));
783:       if (ctx) PetscCall(PetscFree(ctx));
784:       PetscCall(DMSetApplicationContext(first, NULL));
785:     }
786:     for (k = 1; k < nref; k++) {
787:       DM dm = dms[d * nref + k];
788:       if (dm) PetscCall(DMSetApplicationContext(dm, NULL));
789:     }
790:   }

792:   /* destroy DMs */
793:   for (k = 0; k < nd * nref; k++) {
794:     if (dms[k]) PetscCall(DMDestroyShellDMDA(&dms[k]));
795:   }
796:   PetscCall(PetscFree(dms));

798:   /* destroy communicators */
799:   for (k = 0; k < nd; k++) {
800:     if (comms[k] != MPI_COMM_NULL) PetscCall(PetscCommDestroy(&comms[k]));
801:   }
802:   PetscCall(PetscFree(comms));
803:   PetscFunctionReturn(PETSC_SUCCESS);
804: }

806: PetscErrorCode test_basic(void)
807: {
808:   DM           dmF, dmdaF = NULL, dmC = NULL;
809:   Mat          A;
810:   Vec          x, b;
811:   KSP          ksp;
812:   PC           pc;
813:   UserContext *user = NULL;

815:   PetscFunctionBeginUser;
816:   PetscCall(UserContextCreate(PETSC_COMM_WORLD, &user));
817:   PetscCall(HierarchyCreate_Basic(&dmF, &dmC, user));
818:   PetscCall(DMShellGetContext(dmF, &dmdaF));

820:   PetscCall(DMCreateMatrix(dmF, &A));
821:   PetscCall(DMCreateGlobalVector(dmF, &x));
822:   PetscCall(DMCreateGlobalVector(dmF, &b));
823:   PetscCall(ComputeRHS_DMDA(dmdaF, b, user));

825:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
826:   PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix_ShellDA, user));
827:   /*PetscCall(KSPSetOperators(ksp,A,A));*/
828:   PetscCall(KSPSetDM(ksp, dmF));
829:   PetscCall(KSPSetDMActive(ksp, PETSC_TRUE));
830:   PetscCall(KSPSetFromOptions(ksp));
831:   PetscCall(KSPGetPC(ksp, &pc));
832:   PetscCall(PCTelescopeSetUseCoarseDM(pc, PETSC_TRUE));

834:   PetscCall(KSPSolve(ksp, b, x));

836:   if (dmC) PetscCall(DMDestroyShellDMDA(&dmC));
837:   PetscCall(DMDestroyShellDMDA(&dmF));
838:   PetscCall(KSPDestroy(&ksp));
839:   PetscCall(MatDestroy(&A));
840:   PetscCall(VecDestroy(&b));
841:   PetscCall(VecDestroy(&x));
842:   PetscCall(PetscFree(user));
843:   PetscFunctionReturn(PETSC_SUCCESS);
844: }

846: PetscErrorCode test_mg(void)
847: {
848:   DM           dmF, dmdaF = NULL, *dms = NULL;
849:   Mat          A;
850:   Vec          x, b;
851:   KSP          ksp;
852:   PetscInt     k, d, nd, nref;
853:   MPI_Comm    *comms = NULL;
854:   UserContext *user  = NULL;

856:   PetscFunctionBeginUser;
857:   PetscCall(HierarchyCreate(&nd, &nref, &comms, &dms));
858:   dmF = dms[nd * nref - 1];

860:   PetscCall(DMShellGetContext(dmF, &dmdaF));
861:   PetscCall(DMGetApplicationContext(dmF, &user));

863:   PetscCall(DMCreateMatrix(dmF, &A));
864:   PetscCall(DMCreateGlobalVector(dmF, &x));
865:   PetscCall(DMCreateGlobalVector(dmF, &b));
866:   PetscCall(ComputeRHS_DMDA(dmdaF, b, user));

868:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
869:   PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix_ShellDA, user));
870:   /*PetscCall(KSPSetOperators(ksp,A,A));*/
871:   PetscCall(KSPSetDM(ksp, dmF));
872:   PetscCall(KSPSetDMActive(ksp, PETSC_TRUE));
873:   PetscCall(KSPSetFromOptions(ksp));

875:   PetscCall(KSPSolve(ksp, b, x));

877:   for (d = 0; d < nd; d++) {
878:     DM first = dms[d * nref + 0];

880:     if (first) {
881:       UserContext *ctx = NULL;

883:       PetscCall(DMGetApplicationContext(first, &ctx));
884:       if (ctx) PetscCall(PetscFree(ctx));
885:       PetscCall(DMSetApplicationContext(first, NULL));
886:     }
887:     for (k = 1; k < nref; k++) {
888:       DM dm = dms[d * nref + k];
889:       if (dm) PetscCall(DMSetApplicationContext(dm, NULL));
890:     }
891:   }

893:   for (k = 0; k < nd * nref; k++) {
894:     if (dms[k]) PetscCall(DMDestroyShellDMDA(&dms[k]));
895:   }
896:   PetscCall(PetscFree(dms));

898:   PetscCall(KSPDestroy(&ksp));
899:   PetscCall(MatDestroy(&A));
900:   PetscCall(VecDestroy(&b));
901:   PetscCall(VecDestroy(&x));

903:   for (k = 0; k < nd; k++) {
904:     if (comms[k] != MPI_COMM_NULL) PetscCall(PetscCommDestroy(&comms[k]));
905:   }
906:   PetscCall(PetscFree(comms));
907:   PetscFunctionReturn(PETSC_SUCCESS);
908: }

910: int main(int argc, char **argv)
911: {
912:   PetscInt test_id = 0;

914:   PetscFunctionBeginUser;
915:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
916:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-tid", &test_id, NULL));
917:   switch (test_id) {
918:   case 0:
919:     PetscCall(test_basic());
920:     break;
921:   case 1:
922:     PetscCall(test_hierarchy());
923:     break;
924:   case 2:
925:     PetscCall(test_mg());
926:     break;
927:   default:
928:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "-tid must be 0,1,2");
929:   }
930:   PetscCall(PetscFinalize());
931:   return 0;
932: }

934: PetscErrorCode ComputeRHS_DMDA(DM da, Vec b, void *ctx)
935: {
936:   UserContext  *user = (UserContext *)ctx;
937:   PetscInt      i, j, mx, my, xm, ym, xs, ys;
938:   PetscScalar   Hx, Hy;
939:   PetscScalar **array;
940:   PetscBool     isda = PETSC_FALSE;

942:   PetscFunctionBeginUser;
943:   PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isda));
944:   PetscCheck(isda, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "DM provided must be a DMDA");
945:   PetscCall(DMDAGetInfo(da, NULL, &mx, &my, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
946:   Hx = 1.0 / (PetscReal)(mx - 1);
947:   Hy = 1.0 / (PetscReal)(my - 1);
948:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
949:   PetscCall(DMDAVecGetArray(da, b, &array));
950:   for (j = ys; j < ys + ym; j++) {
951:     for (i = xs; i < xs + xm; i++) array[j][i] = PetscExpScalar(-((PetscReal)i * Hx) * ((PetscReal)i * Hx) / user->nu) * PetscExpScalar(-((PetscReal)j * Hy) * ((PetscReal)j * Hy) / user->nu) * Hx * Hy;
952:   }
953:   PetscCall(DMDAVecRestoreArray(da, b, &array));
954:   PetscCall(VecAssemblyBegin(b));
955:   PetscCall(VecAssemblyEnd(b));

957:   /* force right hand side to be consistent for singular matrix */
958:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
959:   if (user->bcType == NEUMANN) {
960:     MatNullSpace nullspace;

962:     PetscCall(MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace));
963:     PetscCall(MatNullSpaceRemove(nullspace, b));
964:     PetscCall(MatNullSpaceDestroy(&nullspace));
965:   }
966:   PetscFunctionReturn(PETSC_SUCCESS);
967: }

969: PetscErrorCode ComputeRho(PetscInt i, PetscInt j, PetscInt mx, PetscInt my, PetscReal centerRho, PetscReal *rho)
970: {
971:   PetscFunctionBeginUser;
972:   if ((i > mx / 3.0) && (i < 2.0 * mx / 3.0) && (j > my / 3.0) && (j < 2.0 * my / 3.0)) {
973:     *rho = centerRho;
974:   } else {
975:     *rho = 1.0;
976:   }
977:   PetscFunctionReturn(PETSC_SUCCESS);
978: }

980: PetscErrorCode ComputeMatrix_DMDA(DM da, Mat J, Mat jac, void *ctx)
981: {
982:   UserContext *user = (UserContext *)ctx;
983:   PetscReal    centerRho;
984:   PetscInt     i, j, mx, my, xm, ym, xs, ys;
985:   PetscScalar  v[5];
986:   PetscReal    Hx, Hy, HydHx, HxdHy, rho;
987:   MatStencil   row, col[5];
988:   PetscBool    isda = PETSC_FALSE;

990:   PetscFunctionBeginUser;
991:   PetscCall(PetscObjectTypeCompare((PetscObject)da, DMDA, &isda));
992:   PetscCheck(isda, PetscObjectComm((PetscObject)da), PETSC_ERR_USER, "DM provided must be a DMDA");
993:   PetscCall(MatZeroEntries(jac));
994:   centerRho = user->rho;
995:   PetscCall(DMDAGetInfo(da, NULL, &mx, &my, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
996:   Hx    = 1.0 / (PetscReal)(mx - 1);
997:   Hy    = 1.0 / (PetscReal)(my - 1);
998:   HxdHy = Hx / Hy;
999:   HydHx = Hy / Hx;
1000:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));
1001:   for (j = ys; j < ys + ym; j++) {
1002:     for (i = xs; i < xs + xm; i++) {
1003:       row.i = i;
1004:       row.j = j;
1005:       PetscCall(ComputeRho(i, j, mx, my, centerRho, &rho));
1006:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
1007:         if (user->bcType == DIRICHLET) {
1008:           v[0] = 2.0 * rho * (HxdHy + HydHx);
1009:           PetscCall(MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES));
1010:         } else if (user->bcType == NEUMANN) {
1011:           PetscInt numx = 0, numy = 0, num = 0;
1012:           if (j != 0) {
1013:             v[num]     = -rho * HxdHy;
1014:             col[num].i = i;
1015:             col[num].j = j - 1;
1016:             numy++;
1017:             num++;
1018:           }
1019:           if (i != 0) {
1020:             v[num]     = -rho * HydHx;
1021:             col[num].i = i - 1;
1022:             col[num].j = j;
1023:             numx++;
1024:             num++;
1025:           }
1026:           if (i != mx - 1) {
1027:             v[num]     = -rho * HydHx;
1028:             col[num].i = i + 1;
1029:             col[num].j = j;
1030:             numx++;
1031:             num++;
1032:           }
1033:           if (j != my - 1) {
1034:             v[num]     = -rho * HxdHy;
1035:             col[num].i = i;
1036:             col[num].j = j + 1;
1037:             numy++;
1038:             num++;
1039:           }
1040:           v[num]     = numx * rho * HydHx + numy * rho * HxdHy;
1041:           col[num].i = i;
1042:           col[num].j = j;
1043:           num++;
1044:           PetscCall(MatSetValuesStencil(jac, 1, &row, num, col, v, INSERT_VALUES));
1045:         }
1046:       } else {
1047:         v[0]     = -rho * HxdHy;
1048:         col[0].i = i;
1049:         col[0].j = j - 1;
1050:         v[1]     = -rho * HydHx;
1051:         col[1].i = i - 1;
1052:         col[1].j = j;
1053:         v[2]     = 2.0 * rho * (HxdHy + HydHx);
1054:         col[2].i = i;
1055:         col[2].j = j;
1056:         v[3]     = -rho * HydHx;
1057:         col[3].i = i + 1;
1058:         col[3].j = j;
1059:         v[4]     = -rho * HxdHy;
1060:         col[4].i = i;
1061:         col[4].j = j + 1;
1062:         PetscCall(MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES));
1063:       }
1064:     }
1065:   }
1066:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
1067:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
1068:   PetscFunctionReturn(PETSC_SUCCESS);
1069: }

1071: PetscErrorCode ComputeMatrix_ShellDA(KSP ksp, Mat J, Mat jac, void *ctx)
1072: {
1073:   DM dm, da;
1074:   PetscFunctionBeginUser;
1075:   PetscCall(KSPGetDM(ksp, &dm));
1076:   PetscCall(DMShellGetContext(dm, &da));
1077:   PetscCall(ComputeMatrix_DMDA(da, J, jac, ctx));
1078:   PetscFunctionReturn(PETSC_SUCCESS);
1079: }

1081: /*TEST

1083:   test:
1084:     suffix: basic_dirichlet
1085:     nsize: 4
1086:     args: -tid 0 -ksp_monitor_short -pc_type telescope -telescope_ksp_max_it 100000 -telescope_pc_type lu -telescope_ksp_type fgmres -telescope_ksp_monitor_short -ksp_type gcr

1088:   test:
1089:     suffix: basic_neumann
1090:     nsize: 4
1091:     requires: !single
1092:     args: -tid 0 -ksp_monitor_short -pc_type telescope -telescope_ksp_max_it 100000 -telescope_pc_type jacobi -telescope_ksp_type fgmres -telescope_ksp_monitor_short -ksp_type gcr -bc_type neumann

1094:   test:
1095:     suffix: mg_2lv_2mg
1096:     nsize: 6
1097:     args: -tid 2 -ksp_monitor_short -ksp_type gcr -pc_type mg -pc_mg_levels 3 -level_nrefs 3 -ndecomps 2 -mg_coarse_pc_type telescope -mg_coarse_pc_telescope_use_coarse_dm -mg_coarse_telescope_pc_type mg -mg_coarse_telescope_pc_mg_levels 3 -level_comm_red_factor 6 -mg_coarse_telescope_mg_coarse_pc_type lu

1099:   test:
1100:     suffix: mg_3lv_2mg
1101:     nsize: 4
1102:     args: -tid 2 -ksp_monitor_short -ksp_type gcr -pc_type mg -pc_mg_levels 3 -level_nrefs 3 -ndecomps 3 -mg_coarse_pc_type telescope -mg_coarse_pc_telescope_use_coarse_dm -mg_coarse_telescope_pc_type mg -mg_coarse_telescope_pc_mg_levels 3 -mg_coarse_telescope_mg_coarse_pc_type telescope -mg_coarse_telescope_mg_coarse_pc_telescope_use_coarse_dm -mg_coarse_telescope_mg_coarse_telescope_pc_type lu -m 5

1104:   test:
1105:     suffix: mg_3lv_2mg_customcommsize
1106:     nsize: 12
1107:     args: -tid 2 -ksp_monitor_short -ksp_type gcr -pc_type mg -pc_mg_levels 3 -level_nrefs 3 -ndecomps 3 -mg_coarse_pc_type telescope -mg_coarse_pc_telescope_use_coarse_dm  -m 5 -level_comm_red_factor 2,6 -mg_coarse_telescope_pc_type mg -mg_coarse_telescope_pc_mg_levels 3 -mg_coarse_telescope_mg_coarse_pc_type telescope -mg_coarse_telescope_mg_coarse_pc_telescope_use_coarse_dm -mg_coarse_telescope_mg_coarse_telescope_pc_type lu

1109:  TEST*/