Actual source code: ex73.c
2: /*
3: This example was derived from src/ksp/ksp/tutorials ex29.c
5: Inhomogeneous Laplacian in 2D. Modeled by the partial differential equation
7: -div \rho grad u = f, 0 < x,y < 1,
9: with forcing function
11: f = e^{-x^2/\nu} e^{-y^2/\nu}
13: with Dirichlet boundary conditions
15: u = f(x,y) for x = 0, x = 1, y = 0, y = 1
17: or pure Neumman boundary conditions.
18: */
20: static char help[] = "Solves 2D inhomogeneous Laplacian. Demonstates using PCTelescopeSetCoarseDM functionality of PCTelescope via a DMShell\n\n";
22: #include <petscdm.h>
23: #include <petscdmda.h>
24: #include <petscdmshell.h>
25: #include <petscksp.h>
27: PetscErrorCode ComputeMatrix_ShellDA(KSP,Mat,Mat,void*);
28: PetscErrorCode ComputeMatrix_DMDA(DM,Mat,Mat,void*);
29: PetscErrorCode ComputeRHS_DMDA(DM,Vec,void*);
30: PetscErrorCode DMShellCreate_ShellDA(DM,DM*);
31: PetscErrorCode DMFieldScatter_ShellDA(DM,Vec,ScatterMode,DM,Vec);
32: PetscErrorCode DMStateScatter_ShellDA(DM,ScatterMode,DM);
34: typedef enum {DIRICHLET, NEUMANN} BCType;
36: typedef struct {
37: PetscReal rho;
38: PetscReal nu;
39: BCType bcType;
40: MPI_Comm comm;
41: } UserContext;
43: PetscErrorCode UserContextCreate(MPI_Comm comm,UserContext **ctx)
44: {
45: UserContext *user;
46: const char *bcTypes[2] = {"dirichlet","neumann"};
47: PetscInt bc;
51: PetscCalloc1(1,&user);
52: user->comm = comm;
53: PetscOptionsBegin(comm, "", "Options for the inhomogeneous Poisson equation", "DMqq");
54: user->rho = 1.0;
55: PetscOptionsReal("-rho", "The conductivity", "ex29.c", user->rho, &user->rho, NULL);
56: user->nu = 0.1;
57: PetscOptionsReal("-nu", "The width of the Gaussian source", "ex29.c", user->nu, &user->nu, NULL);
58: bc = (PetscInt)DIRICHLET;
59: PetscOptionsEList("-bc_type","Type of boundary condition","ex29.c",bcTypes,2,bcTypes[0],&bc,NULL);
60: user->bcType = (BCType)bc;
61: PetscOptionsEnd();
62: *ctx = user;
63: return 0;
64: }
66: PetscErrorCode CommCoarsen(MPI_Comm comm,PetscInt number,PetscSubcomm *p)
67: {
68: PetscSubcomm psubcomm;
70: PetscSubcommCreate(comm,&psubcomm);
71: PetscSubcommSetNumber(psubcomm,number);
72: PetscSubcommSetType(psubcomm,PETSC_SUBCOMM_INTERLACED);
73: *p = psubcomm;
74: return 0;
75: }
77: PetscErrorCode CommHierarchyCreate(MPI_Comm comm,PetscInt n,PetscInt number[],PetscSubcomm pscommlist[])
78: {
79: PetscInt k;
80: PetscBool view_hierarchy = PETSC_FALSE;
83: for (k=0; k<n; k++) {
84: pscommlist[k] = NULL;
85: }
87: if (n < 1) return 0;
89: CommCoarsen(comm,number[n-1],&pscommlist[n-1]);
90: for (k=n-2; k>=0; k--) {
91: MPI_Comm comm_k = PetscSubcommChild(pscommlist[k+1]);
92: if (pscommlist[k+1]->color == 0) {
93: CommCoarsen(comm_k,number[k],&pscommlist[k]);
94: }
95: }
97: PetscOptionsGetBool(NULL,NULL,"-view_hierarchy",&view_hierarchy,NULL);
98: if (view_hierarchy) {
99: PetscMPIInt size;
101: MPI_Comm_size(comm,&size);
102: PetscPrintf(comm,"level[%D] size %d\n",n,(int)size);
103: for (k=n-1; k>=0; k--) {
104: if (pscommlist[k]) {
105: MPI_Comm comm_k = PetscSubcommChild(pscommlist[k]);
107: if (pscommlist[k]->color == 0) {
108: MPI_Comm_size(comm_k,&size);
109: PetscPrintf(comm_k,"level[%D] size %d\n",k,(int)size);
110: }
111: }
112: }
113: }
114: return 0;
115: }
117: /* taken from src/ksp/pc/impls/telescope/telescope_dmda.c */
118: static PetscErrorCode _DMDADetermineRankFromGlobalIJ_2d(PetscInt i,PetscInt j,PetscInt Mp,PetscInt Np,
119: PetscInt start_i[],PetscInt start_j[],
120: PetscInt span_i[],PetscInt span_j[],
121: PetscMPIInt *_pi,PetscMPIInt *_pj,PetscMPIInt *rank_re)
122: {
123: PetscInt pi,pj,n;
126: *rank_re = -1;
127: pi = pj = -1;
128: if (_pi) {
129: for (n=0; n<Mp; n++) {
130: if ((i >= start_i[n]) && (i < start_i[n]+span_i[n])) {
131: pi = n;
132: break;
133: }
134: }
136: *_pi = (PetscMPIInt)pi;
137: }
139: if (_pj) {
140: for (n=0; n<Np; n++) {
141: if ((j >= start_j[n]) && (j < start_j[n]+span_j[n])) {
142: pj = n;
143: break;
144: }
145: }
147: *_pj = (PetscMPIInt)pj;
148: }
150: *rank_re = (PetscMPIInt)(pi + pj * Mp);
151: return 0;
152: }
154: /* taken from src/ksp/pc/impls/telescope/telescope_dmda.c */
155: static PetscErrorCode _DMDADetermineGlobalS0_2d(PetscMPIInt rank_re,PetscInt Mp_re,PetscInt Np_re,
156: PetscInt range_i_re[],PetscInt range_j_re[],PetscInt *s0)
157: {
158: PetscInt i,j,start_IJ = 0;
159: PetscMPIInt rank_ij;
162: *s0 = -1;
163: for (j=0; j<Np_re; j++) {
164: for (i=0; i<Mp_re; i++) {
165: rank_ij = (PetscMPIInt)(i + j*Mp_re);
166: if (rank_ij < rank_re) {
167: start_IJ += range_i_re[i]*range_j_re[j];
168: }
169: }
170: }
171: *s0 = start_IJ;
172: return 0;
173: }
175: /* adapted from src/ksp/pc/impls/telescope/telescope_dmda.c */
176: static PetscErrorCode DMDACreatePermutation_2d(DM dmrepart,DM dmf,Mat *mat)
177: {
179: PetscInt k,sum,Mp_re = 0,Np_re = 0;
180: PetscInt nx,ny,sr,er,Mr,ndof;
181: PetscInt i,j,location,startI[2],endI[2],lenI[2];
182: const PetscInt *_range_i_re = NULL,*_range_j_re = NULL;
183: PetscInt *range_i_re,*range_j_re;
184: PetscInt *start_i_re,*start_j_re;
185: MPI_Comm comm;
186: Vec V;
187: Mat Pscalar;
190: PetscInfo(dmf,"setting up the permutation matrix (DMDA-2D)\n");
191: PetscObjectGetComm((PetscObject)dmf,&comm);
193: _range_i_re = _range_j_re = NULL;
194: /* Create DMDA on the child communicator */
195: if (dmrepart) {
196: DMDAGetInfo(dmrepart,NULL,NULL,NULL,NULL,&Mp_re,&Np_re,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
197: DMDAGetOwnershipRanges(dmrepart,&_range_i_re,&_range_j_re,NULL);
198: }
200: /* note - assume rank 0 always participates */
201: MPI_Bcast(&Mp_re,1,MPIU_INT,0,comm);
202: MPI_Bcast(&Np_re,1,MPIU_INT,0,comm);
204: PetscCalloc1(Mp_re,&range_i_re);
205: PetscCalloc1(Np_re,&range_j_re);
207: if (_range_i_re) PetscArraycpy(range_i_re,_range_i_re,Mp_re);
208: if (_range_j_re) PetscArraycpy(range_j_re,_range_j_re,Np_re);
210: MPI_Bcast(range_i_re,Mp_re,MPIU_INT,0,comm);
211: MPI_Bcast(range_j_re,Np_re,MPIU_INT,0,comm);
213: PetscMalloc1(Mp_re,&start_i_re);
214: PetscMalloc1(Np_re,&start_j_re);
216: sum = 0;
217: for (k=0; k<Mp_re; k++) {
218: start_i_re[k] = sum;
219: sum += range_i_re[k];
220: }
222: sum = 0;
223: for (k=0; k<Np_re; k++) {
224: start_j_re[k] = sum;
225: sum += range_j_re[k];
226: }
228: /* Create permutation */
229: DMDAGetInfo(dmf,NULL,&nx,&ny,NULL,NULL,NULL,NULL,&ndof,NULL,NULL,NULL,NULL,NULL);
230: DMGetGlobalVector(dmf,&V);
231: VecGetSize(V,&Mr);
232: VecGetOwnershipRange(V,&sr,&er);
233: DMRestoreGlobalVector(dmf,&V);
234: sr = sr / ndof;
235: er = er / ndof;
236: Mr = Mr / ndof;
238: MatCreate(comm,&Pscalar);
239: MatSetSizes(Pscalar,(er-sr),(er-sr),Mr,Mr);
240: MatSetType(Pscalar,MATAIJ);
241: MatSeqAIJSetPreallocation(Pscalar,1,NULL);
242: MatMPIAIJSetPreallocation(Pscalar,1,NULL,1,NULL);
244: DMDAGetCorners(dmf,NULL,NULL,NULL,&lenI[0],&lenI[1],NULL);
245: DMDAGetCorners(dmf,&startI[0],&startI[1],NULL,&endI[0],&endI[1],NULL);
246: endI[0] += startI[0];
247: endI[1] += startI[1];
249: for (j=startI[1]; j<endI[1]; j++) {
250: for (i=startI[0]; i<endI[0]; i++) {
251: PetscMPIInt rank_ijk_re,rank_reI[] = {0,0};
252: PetscInt s0_re;
253: PetscInt ii,jj,local_ijk_re,mapped_ijk;
254: PetscInt lenI_re[] = {0,0};
256: location = (i - startI[0]) + (j - startI[1])*lenI[0];
257: _DMDADetermineRankFromGlobalIJ_2d(i,j,Mp_re,Np_re,
258: start_i_re,start_j_re,
259: range_i_re,range_j_re,
260: &rank_reI[0],&rank_reI[1],&rank_ijk_re);
262: _DMDADetermineGlobalS0_2d(rank_ijk_re,Mp_re,Np_re,range_i_re,range_j_re,&s0_re);
264: ii = i - start_i_re[ rank_reI[0] ];
266: jj = j - start_j_re[ rank_reI[1] ];
269: lenI_re[0] = range_i_re[ rank_reI[0] ];
270: lenI_re[1] = range_j_re[ rank_reI[1] ];
271: local_ijk_re = ii + jj * lenI_re[0];
272: mapped_ijk = s0_re + local_ijk_re;
273: MatSetValue(Pscalar,sr+location,mapped_ijk,1.0,INSERT_VALUES);
274: }
275: }
276: MatAssemblyBegin(Pscalar,MAT_FINAL_ASSEMBLY);
277: MatAssemblyEnd(Pscalar,MAT_FINAL_ASSEMBLY);
279: *mat = Pscalar;
281: PetscFree(range_i_re);
282: PetscFree(range_j_re);
283: PetscFree(start_i_re);
284: PetscFree(start_j_re);
285: return 0;
286: }
288: /* adapted from src/ksp/pc/impls/telescope/telescope_dmda.c */
289: static PetscErrorCode PCTelescopeSetUp_dmda_scatters(DM dmf,DM dmc)
290: {
291: Vec xred,yred,xtmp,x,xp;
292: VecScatter scatter;
293: IS isin;
294: PetscInt m,bs,st,ed;
295: MPI_Comm comm;
296: VecType vectype;
299: PetscObjectGetComm((PetscObject)dmf,&comm);
300: DMCreateGlobalVector(dmf,&x);
301: VecGetBlockSize(x,&bs);
302: VecGetType(x,&vectype);
304: /* cannot use VecDuplicate as xp is already composed with dmf */
305: /*VecDuplicate(x,&xp);*/
306: {
307: PetscInt m,M;
309: VecGetSize(x,&M);
310: VecGetLocalSize(x,&m);
311: VecCreate(comm,&xp);
312: VecSetSizes(xp,m,M);
313: VecSetBlockSize(xp,bs);
314: VecSetType(xp,vectype);
315: }
317: m = 0;
318: xred = NULL;
319: yred = NULL;
320: if (dmc) {
321: DMCreateGlobalVector(dmc,&xred);
322: VecDuplicate(xred,&yred);
323: VecGetOwnershipRange(xred,&st,&ed);
324: ISCreateStride(comm,ed-st,st,1,&isin);
325: VecGetLocalSize(xred,&m);
326: } else {
327: VecGetOwnershipRange(x,&st,&ed);
328: ISCreateStride(comm,0,st,1,&isin);
329: }
330: ISSetBlockSize(isin,bs);
331: VecCreate(comm,&xtmp);
332: VecSetSizes(xtmp,m,PETSC_DECIDE);
333: VecSetBlockSize(xtmp,bs);
334: VecSetType(xtmp,vectype);
335: VecScatterCreate(x,isin,xtmp,NULL,&scatter);
337: PetscObjectCompose((PetscObject)dmf,"isin",(PetscObject)isin);
338: PetscObjectCompose((PetscObject)dmf,"scatter",(PetscObject)scatter);
339: PetscObjectCompose((PetscObject)dmf,"xtmp",(PetscObject)xtmp);
340: PetscObjectCompose((PetscObject)dmf,"xp",(PetscObject)xp);
342: VecDestroy(&xred);
343: VecDestroy(&yred);
344: VecDestroy(&x);
345: return 0;
346: }
348: PetscErrorCode DMCreateMatrix_ShellDA(DM dm,Mat *A)
349: {
350: DM da;
351: MPI_Comm comm;
352: PetscMPIInt size;
353: UserContext *ctx = NULL;
354: PetscInt M,N;
357: DMShellGetContext(dm,&da);
358: PetscObjectGetComm((PetscObject)da,&comm);
359: MPI_Comm_size(comm,&size);
360: DMCreateMatrix(da,A);
361: MatGetSize(*A,&M,&N);
362: PetscPrintf(comm,"[size %D] DMCreateMatrix_ShellDA (%D x %D)\n",(PetscInt)size,M,N);
364: DMGetApplicationContext(dm,&ctx);
365: if (ctx->bcType == NEUMANN) {
366: MatNullSpace nullspace = NULL;
367: PetscPrintf(comm,"[size %D] DMCreateMatrix_ShellDA: using neumann bcs\n",(PetscInt)size);
369: MatGetNullSpace(*A,&nullspace);
370: if (!nullspace) {
371: PetscPrintf(comm,"[size %D] DMCreateMatrix_ShellDA: operator does not have nullspace - attaching\n",(PetscInt)size);
372: MatNullSpaceCreate(comm,PETSC_TRUE,0,0,&nullspace);
373: MatSetNullSpace(*A,nullspace);
374: MatNullSpaceDestroy(&nullspace);
375: } else {
376: PetscPrintf(comm,"[size %D] DMCreateMatrix_ShellDA: operator already has a nullspace\n",(PetscInt)size);
377: }
378: }
379: return 0;
380: }
382: PetscErrorCode DMCreateGlobalVector_ShellDA(DM dm,Vec *x)
383: {
384: DM da;
386: DMShellGetContext(dm,&da);
387: DMCreateGlobalVector(da,x);
388: VecSetDM(*x,dm);
389: return 0;
390: }
392: PetscErrorCode DMCreateLocalVector_ShellDA(DM dm,Vec *x)
393: {
394: DM da;
396: DMShellGetContext(dm,&da);
397: DMCreateLocalVector(da,x);
398: VecSetDM(*x,dm);
399: return 0;
400: }
402: PetscErrorCode DMCoarsen_ShellDA(DM dm,MPI_Comm comm,DM *dmc)
403: {
405: *dmc = NULL;
406: DMGetCoarseDM(dm,dmc);
407: if (!*dmc) {
408: SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"The coarse DM should never be NULL. The DM hierarchy should have already been defined");
409: } else {
410: PetscObjectReference((PetscObject)(*dmc));
411: }
412: return 0;
413: }
415: PetscErrorCode DMCreateInterpolation_ShellDA(DM dm1,DM dm2,Mat *mat,Vec *vec)
416: {
417: DM da1,da2;
419: DMShellGetContext(dm1,&da1);
420: DMShellGetContext(dm2,&da2);
421: DMCreateInterpolation(da1,da2,mat,vec);
422: return 0;
423: }
425: PetscErrorCode DMShellDASetUp_TelescopeDMScatter(DM dmf_shell,DM dmc_shell)
426: {
427: Mat P = NULL;
428: DM dmf = NULL,dmc = NULL;
431: DMShellGetContext(dmf_shell,&dmf);
432: if (dmc_shell) {
433: DMShellGetContext(dmc_shell,&dmc);
434: }
435: DMDACreatePermutation_2d(dmc,dmf,&P);
436: PetscObjectCompose((PetscObject)dmf,"P",(PetscObject)P);
437: PCTelescopeSetUp_dmda_scatters(dmf,dmc);
438: PetscObjectComposeFunction((PetscObject)dmf_shell,"PCTelescopeFieldScatter",DMFieldScatter_ShellDA);
439: PetscObjectComposeFunction((PetscObject)dmf_shell,"PCTelescopeStateScatter",DMStateScatter_ShellDA);
440: return 0;
441: }
443: PetscErrorCode DMShellDAFieldScatter_Forward(DM dmf,Vec x,DM dmc,Vec xc)
444: {
445: Mat P = NULL;
446: Vec xp = NULL,xtmp = NULL;
447: VecScatter scatter = NULL;
448: const PetscScalar *x_array;
449: PetscInt i,st,ed;
452: PetscObjectQuery((PetscObject)dmf,"P",(PetscObject*)&P);
453: PetscObjectQuery((PetscObject)dmf,"xp",(PetscObject*)&xp);
454: PetscObjectQuery((PetscObject)dmf,"scatter",(PetscObject*)&scatter);
455: PetscObjectQuery((PetscObject)dmf,"xtmp",(PetscObject*)&xtmp);
461: MatMultTranspose(P,x,xp);
463: /* pull in vector x->xtmp */
464: VecScatterBegin(scatter,xp,xtmp,INSERT_VALUES,SCATTER_FORWARD);
465: VecScatterEnd(scatter,xp,xtmp,INSERT_VALUES,SCATTER_FORWARD);
467: /* copy vector entries into xred */
468: VecGetArrayRead(xtmp,&x_array);
469: if (xc) {
470: PetscScalar *LA_xred;
471: VecGetOwnershipRange(xc,&st,&ed);
473: VecGetArray(xc,&LA_xred);
474: for (i=0; i<ed-st; i++) {
475: LA_xred[i] = x_array[i];
476: }
477: VecRestoreArray(xc,&LA_xred);
478: }
479: VecRestoreArrayRead(xtmp,&x_array);
480: return 0;
481: }
483: PetscErrorCode DMShellDAFieldScatter_Reverse(DM dmf,Vec y,DM dmc,Vec yc)
484: {
485: Mat P = NULL;
486: Vec xp = NULL,xtmp = NULL;
487: VecScatter scatter = NULL;
488: PetscScalar *array;
489: PetscInt i,st,ed;
492: PetscObjectQuery((PetscObject)dmf,"P",(PetscObject*)&P);
493: PetscObjectQuery((PetscObject)dmf,"xp",(PetscObject*)&xp);
494: PetscObjectQuery((PetscObject)dmf,"scatter",(PetscObject*)&scatter);
495: PetscObjectQuery((PetscObject)dmf,"xtmp",(PetscObject*)&xtmp);
502: /* return vector */
503: VecGetArray(xtmp,&array);
504: if (yc) {
505: const PetscScalar *LA_yred;
506: VecGetOwnershipRange(yc,&st,&ed);
507: VecGetArrayRead(yc,&LA_yred);
508: for (i=0; i<ed-st; i++) {
509: array[i] = LA_yred[i];
510: }
511: VecRestoreArrayRead(yc,&LA_yred);
512: }
513: VecRestoreArray(xtmp,&array);
514: VecScatterBegin(scatter,xtmp,xp,INSERT_VALUES,SCATTER_REVERSE);
515: VecScatterEnd(scatter,xtmp,xp,INSERT_VALUES,SCATTER_REVERSE);
516: MatMult(P,xp,y);
517: return 0;
518: }
520: PetscErrorCode DMFieldScatter_ShellDA(DM dmf_shell,Vec x,ScatterMode mode,DM dmc_shell,Vec xc)
521: {
522: DM dmf = NULL,dmc = NULL;
525: DMShellGetContext(dmf_shell,&dmf);
526: if (dmc_shell) {
527: DMShellGetContext(dmc_shell,&dmc);
528: }
529: if (mode == SCATTER_FORWARD) {
530: DMShellDAFieldScatter_Forward(dmf,x,dmc,xc);
531: } else if (mode == SCATTER_REVERSE) {
532: DMShellDAFieldScatter_Reverse(dmf,x,dmc,xc);
533: } else SETERRQ(PetscObjectComm((PetscObject)dmf_shell),PETSC_ERR_SUP,"Only mode = SCATTER_FORWARD, SCATTER_REVERSE supported");
534: return 0;
535: }
537: PetscErrorCode DMStateScatter_ShellDA(DM dmf_shell,ScatterMode mode,DM dmc_shell)
538: {
539: PetscMPIInt size_f = 0,size_c = 0;
541: MPI_Comm_size(PetscObjectComm((PetscObject)dmf_shell),&size_f);
542: if (dmc_shell) {
543: MPI_Comm_size(PetscObjectComm((PetscObject)dmc_shell),&size_c);
544: }
545: if (mode == SCATTER_FORWARD) {
546: PetscPrintf(PetscObjectComm((PetscObject)dmf_shell),"User supplied state scatter (fine [size %d]-> coarse [size %d])\n",(int)size_f,(int)size_c);
547: } else if (mode == SCATTER_REVERSE) {
548: } else SETERRQ(PetscObjectComm((PetscObject)dmf_shell),PETSC_ERR_SUP,"Only mode = SCATTER_FORWARD, SCATTER_REVERSE supported");
549: return 0;
550: }
552: PetscErrorCode DMShellCreate_ShellDA(DM da,DM *dms)
553: {
555: if (da) {
556: DMShellCreate(PetscObjectComm((PetscObject)da),dms);
557: DMShellSetContext(*dms,da);
558: DMShellSetCreateGlobalVector(*dms,DMCreateGlobalVector_ShellDA);
559: DMShellSetCreateLocalVector(*dms,DMCreateLocalVector_ShellDA);
560: DMShellSetCreateMatrix(*dms,DMCreateMatrix_ShellDA);
561: DMShellSetCoarsen(*dms,DMCoarsen_ShellDA);
562: DMShellSetCreateInterpolation(*dms,DMCreateInterpolation_ShellDA);
563: } else {
564: *dms = NULL;
565: }
566: return 0;
567: }
569: PetscErrorCode DMDestroyShellDMDA(DM *_dm)
570: {
571: DM dm,da = NULL;
574: if (!_dm) return 0;
575: dm = *_dm;
576: if (!dm) return 0;
578: DMShellGetContext(dm,&da);
579: if (da) {
580: Vec vec;
581: VecScatter scatter = NULL;
582: IS is = NULL;
583: Mat P = NULL;
585: PetscObjectQuery((PetscObject)da,"P",(PetscObject*)&P);
586: MatDestroy(&P);
588: vec = NULL;
589: PetscObjectQuery((PetscObject)da,"xp",(PetscObject*)&vec);
590: VecDestroy(&vec);
592: PetscObjectQuery((PetscObject)da,"scatter",(PetscObject*)&scatter);
593: VecScatterDestroy(&scatter);
595: vec = NULL;
596: PetscObjectQuery((PetscObject)da,"xtmp",(PetscObject*)&vec);
597: VecDestroy(&vec);
599: PetscObjectQuery((PetscObject)da,"isin",(PetscObject*)&is);
600: ISDestroy(&is);
602: DMDestroy(&da);
603: }
604: DMDestroy(&dm);
605: *_dm = NULL;
606: return 0;
607: }
609: PetscErrorCode HierarchyCreate_Basic(DM *dm_f,DM *dm_c,UserContext *ctx)
610: {
611: DM dm,dmc,dm_shell,dmc_shell;
612: PetscMPIInt rank;
615: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
616: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,17,17,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&dm);
617: DMSetFromOptions(dm);
618: DMSetUp(dm);
619: DMDASetUniformCoordinates(dm,0,1,0,1,0,0);
620: DMDASetFieldName(dm,0,"Pressure");
621: DMShellCreate_ShellDA(dm,&dm_shell);
622: DMSetApplicationContext(dm_shell,ctx);
624: dmc = NULL;
625: dmc_shell = NULL;
626: if (rank == 0) {
627: DMDACreate2d(PETSC_COMM_SELF,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,17,17,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&dmc);
628: DMSetFromOptions(dmc);
629: DMSetUp(dmc);
630: DMDASetUniformCoordinates(dmc,0,1,0,1,0,0);
631: DMDASetFieldName(dmc,0,"Pressure");
632: DMShellCreate_ShellDA(dmc,&dmc_shell);
633: DMSetApplicationContext(dmc_shell,ctx);
634: }
636: DMSetCoarseDM(dm_shell,dmc_shell);
637: DMShellDASetUp_TelescopeDMScatter(dm_shell,dmc_shell);
639: *dm_f = dm_shell;
640: *dm_c = dmc_shell;
641: return 0;
642: }
644: PetscErrorCode HierarchyCreate(PetscInt *_nd,PetscInt *_nref,MPI_Comm **_cl,DM **_dl)
645: {
646: PetscInt d,k,ndecomps,ncoarsen,found,nx;
647: PetscInt levelrefs,*number;
648: PetscSubcomm *pscommlist;
649: MPI_Comm *commlist;
650: DM *dalist,*dmlist;
651: PetscBool set;
654: ndecomps = 1;
655: PetscOptionsGetInt(NULL,NULL,"-ndecomps",&ndecomps,NULL);
656: ncoarsen = ndecomps - 1;
659: levelrefs = 2;
660: PetscOptionsGetInt(NULL,NULL,"-level_nrefs",&levelrefs,NULL);
661: PetscMalloc1(ncoarsen+1,&number);
662: for (k=0; k<ncoarsen+1; k++) {
663: number[k] = 2;
664: }
665: found = ncoarsen;
666: set = PETSC_FALSE;
667: PetscOptionsGetIntArray(NULL,NULL,"-level_comm_red_factor",number,&found,&set);
668: if (set) {
670: }
672: PetscMalloc1(ncoarsen+1,&pscommlist);
673: for (k=0; k<ncoarsen+1; k++) {
674: pscommlist[k] = NULL;
675: }
677: PetscMalloc1(ndecomps,&commlist);
678: for (k=0; k<ndecomps; k++) {
679: commlist[k] = MPI_COMM_NULL;
680: }
681: PetscMalloc1(ndecomps*levelrefs,&dalist);
682: PetscMalloc1(ndecomps*levelrefs,&dmlist);
683: for (k=0; k<ndecomps*levelrefs; k++) {
684: dalist[k] = NULL;
685: dmlist[k] = NULL;
686: }
688: CommHierarchyCreate(PETSC_COMM_WORLD,ncoarsen,number,pscommlist);
689: for (k=0; k<ncoarsen; k++) {
690: if (pscommlist[k]) {
691: MPI_Comm comm_k = PetscSubcommChild(pscommlist[k]);
692: if (pscommlist[k]->color == 0) {
693: PetscCommDuplicate(comm_k,&commlist[k],NULL);
694: }
695: }
696: }
697: PetscCommDuplicate(PETSC_COMM_WORLD,&commlist[ndecomps-1],NULL);
699: for (k=0; k<ncoarsen; k++) {
700: if (pscommlist[k]) {
701: PetscSubcommDestroy(&pscommlist[k]);
702: }
703: }
705: nx = 17;
706: PetscOptionsGetInt(NULL,NULL,"-m",&nx,NULL);
707: for (d=0; d<ndecomps; d++) {
708: DM dmroot = NULL;
709: char name[PETSC_MAX_PATH_LEN];
711: if (commlist[d] != MPI_COMM_NULL) {
712: DMDACreate2d(commlist[d],DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,nx,nx,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&dmroot);
713: DMSetUp(dmroot);
714: DMDASetUniformCoordinates(dmroot,0,1,0,1,0,0);
715: DMDASetFieldName(dmroot,0,"Pressure");
716: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"root-decomp-%D",d);
717: PetscObjectSetName((PetscObject)dmroot,name);
718: /*DMView(dmroot,PETSC_VIEWER_STDOUT_(commlist[d]));*/
719: }
721: dalist[d*levelrefs + 0] = dmroot;
722: for (k=1; k<levelrefs; k++) {
723: DM dmref = NULL;
725: if (commlist[d] != MPI_COMM_NULL) {
726: DMRefine(dalist[d*levelrefs + (k-1)],MPI_COMM_NULL,&dmref);
727: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"ref%D-decomp-%D",k,d);
728: PetscObjectSetName((PetscObject)dmref,name);
729: DMDAGetInfo(dmref,NULL,&nx,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
730: /*DMView(dmref,PETSC_VIEWER_STDOUT_(commlist[d]));*/
731: }
732: dalist[d*levelrefs + k] = dmref;
733: }
734: MPI_Allreduce(MPI_IN_PLACE,&nx,1,MPIU_INT,MPI_MAX,PETSC_COMM_WORLD);
735: }
737: /* create the hierarchy of DMShell's */
738: for (d=0; d<ndecomps; d++) {
739: char name[PETSC_MAX_PATH_LEN];
741: UserContext *ctx = NULL;
742: if (commlist[d] != MPI_COMM_NULL) {
743: UserContextCreate(commlist[d],&ctx);
744: for (k=0; k<levelrefs; k++) {
745: DMShellCreate_ShellDA(dalist[d*levelrefs + k],&dmlist[d*levelrefs + k]);
746: DMSetApplicationContext(dmlist[d*levelrefs + k],ctx);
747: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"level%D-decomp-%D",k,d);
748: PetscObjectSetName((PetscObject)dmlist[d*levelrefs + k],name);
749: }
750: }
751: }
753: /* set all the coarse DMs */
754: for (k=1; k<ndecomps*levelrefs; k++) { /* skip first DM as it doesn't have a coarse representation */
755: DM dmfine = NULL,dmcoarse = NULL;
757: dmfine = dmlist[k];
758: dmcoarse = dmlist[k-1];
759: if (dmfine) {
760: DMSetCoarseDM(dmfine,dmcoarse);
761: }
762: }
764: /* do special setup on the fine DM coupling different decompositions */
765: for (d=1; d<ndecomps; d++) { /* skip first decomposition as it doesn't have a coarse representation */
766: DM dmfine = NULL,dmcoarse = NULL;
768: dmfine = dmlist[d*levelrefs + 0];
769: dmcoarse = dmlist[(d-1)*levelrefs + (levelrefs-1)];
770: if (dmfine) {
771: DMShellDASetUp_TelescopeDMScatter(dmfine,dmcoarse);
772: }
773: }
775: PetscFree(number);
776: for (k=0; k<ncoarsen; k++) {
777: PetscSubcommDestroy(&pscommlist[k]);
778: }
779: PetscFree(pscommlist);
781: if (_nd) {
782: *_nd = ndecomps;
783: }
784: if (_nref) {
785: *_nref = levelrefs;
786: }
787: if (_cl) {
788: *_cl = commlist;
789: } else {
790: for (k=0; k<ndecomps; k++) {
791: if (commlist[k] != MPI_COMM_NULL) {
792: PetscCommDestroy(&commlist[k]);
793: }
794: }
795: PetscFree(commlist);
796: }
797: if (_dl) {
798: *_dl = dmlist;
799: PetscFree(dalist);
800: } else {
801: for (k=0; k<ndecomps*levelrefs; k++) {
802: DMDestroy(&dmlist[k]);
803: DMDestroy(&dalist[k]);
804: }
805: PetscFree(dmlist);
806: PetscFree(dalist);
807: }
808: return 0;
809: }
811: PetscErrorCode test_hierarchy(void)
812: {
813: PetscInt d,k,nd,nref;
814: MPI_Comm *comms;
815: DM *dms;
818: HierarchyCreate(&nd,&nref,&comms,&dms);
820: /* destroy user context */
821: for (d=0; d<nd; d++) {
822: DM first = dms[d*nref+0];
824: if (first) {
825: UserContext *ctx = NULL;
827: DMGetApplicationContext(first,&ctx);
828: if (ctx) PetscFree(ctx);
829: DMSetApplicationContext(first,NULL);
830: }
831: for (k=1; k<nref; k++) {
832: DM dm = dms[d*nref+k];
833: if (dm) {
834: DMSetApplicationContext(dm,NULL);
835: }
836: }
837: }
839: /* destroy DMs */
840: for (k=0; k<nd*nref; k++) {
841: if (dms[k]) {
842: DMDestroyShellDMDA(&dms[k]);
843: }
844: }
845: PetscFree(dms);
847: /* destroy communicators */
848: for (k=0; k<nd; k++) {
849: if (comms[k] != MPI_COMM_NULL) {
850: PetscCommDestroy(&comms[k]);
851: }
852: }
853: PetscFree(comms);
854: return 0;
855: }
857: PetscErrorCode test_basic(void)
858: {
859: DM dmF,dmdaF = NULL,dmC = NULL;
860: Mat A;
861: Vec x,b;
862: KSP ksp;
863: PC pc;
864: UserContext *user = NULL;
867: UserContextCreate(PETSC_COMM_WORLD,&user);
868: HierarchyCreate_Basic(&dmF,&dmC,user);
869: DMShellGetContext(dmF,&dmdaF);
871: DMCreateMatrix(dmF,&A);
872: DMCreateGlobalVector(dmF,&x);
873: DMCreateGlobalVector(dmF,&b);
874: ComputeRHS_DMDA(dmdaF,b,user);
876: KSPCreate(PETSC_COMM_WORLD,&ksp);
877: KSPSetComputeOperators(ksp,ComputeMatrix_ShellDA,user);
878: /*KSPSetOperators(ksp,A,A);*/
879: KSPSetDM(ksp,dmF);
880: KSPSetDMActive(ksp,PETSC_TRUE);
881: KSPSetFromOptions(ksp);
882: KSPGetPC(ksp,&pc);
883: PCTelescopeSetUseCoarseDM(pc,PETSC_TRUE);
885: KSPSolve(ksp,b,x);
887: if (dmC) {
888: DMDestroyShellDMDA(&dmC);
889: }
890: DMDestroyShellDMDA(&dmF);
891: KSPDestroy(&ksp);
892: MatDestroy(&A);
893: VecDestroy(&b);
894: VecDestroy(&x);
895: PetscFree(user);
896: return 0;
897: }
899: PetscErrorCode test_mg(void)
900: {
901: DM dmF,dmdaF = NULL,*dms = NULL;
902: Mat A;
903: Vec x,b;
904: KSP ksp;
905: PetscInt k,d,nd,nref;
906: MPI_Comm *comms = NULL;
907: UserContext *user = NULL;
910: HierarchyCreate(&nd,&nref,&comms,&dms);
911: dmF = dms[nd*nref-1];
913: DMShellGetContext(dmF,&dmdaF);
914: DMGetApplicationContext(dmF,&user);
916: DMCreateMatrix(dmF,&A);
917: DMCreateGlobalVector(dmF,&x);
918: DMCreateGlobalVector(dmF,&b);
919: ComputeRHS_DMDA(dmdaF,b,user);
921: KSPCreate(PETSC_COMM_WORLD,&ksp);
922: KSPSetComputeOperators(ksp,ComputeMatrix_ShellDA,user);
923: /*KSPSetOperators(ksp,A,A);*/
924: KSPSetDM(ksp,dmF);
925: KSPSetDMActive(ksp,PETSC_TRUE);
926: KSPSetFromOptions(ksp);
928: KSPSolve(ksp,b,x);
930: for (d=0; d<nd; d++) {
931: DM first = dms[d*nref+0];
933: if (first) {
934: UserContext *ctx = NULL;
936: DMGetApplicationContext(first,&ctx);
937: if (ctx) PetscFree(ctx);
938: DMSetApplicationContext(first,NULL);
939: }
940: for (k=1; k<nref; k++) {
941: DM dm = dms[d*nref+k];
942: if (dm) {
943: DMSetApplicationContext(dm,NULL);
944: }
945: }
946: }
948: for (k=0; k<nd*nref; k++) {
949: if (dms[k]) {
950: DMDestroyShellDMDA(&dms[k]);
951: }
952: }
953: PetscFree(dms);
955: KSPDestroy(&ksp);
956: MatDestroy(&A);
957: VecDestroy(&b);
958: VecDestroy(&x);
960: for (k=0; k<nd; k++) {
961: if (comms[k] != MPI_COMM_NULL) {
962: PetscCommDestroy(&comms[k]);
963: }
964: }
965: PetscFree(comms);
966: return 0;
967: }
969: int main(int argc,char **argv)
970: {
971: PetscInt test_id = 0;
973: PetscInitialize(&argc,&argv,(char*)0,help);
974: PetscOptionsGetInt(NULL,NULL,"-tid",&test_id,NULL);
975: switch (test_id) {
976: case 0:
977: test_basic();
978: break;
979: case 1:
980: test_hierarchy();
981: break;
982: case 2:
983: test_mg();
984: break;
985: default:
986: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"-tid must be 0,1,2");
987: }
988: PetscFinalize();
989: return 0;
990: }
992: PetscErrorCode ComputeRHS_DMDA(DM da,Vec b,void *ctx)
993: {
994: UserContext *user = (UserContext*)ctx;
995: PetscInt i,j,mx,my,xm,ym,xs,ys;
996: PetscScalar Hx,Hy;
997: PetscScalar **array;
998: PetscBool isda = PETSC_FALSE;
1001: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
1003: DMDAGetInfo(da,NULL,&mx,&my,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1004: Hx = 1.0 / (PetscReal)(mx-1);
1005: Hy = 1.0 / (PetscReal)(my-1);
1006: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
1007: DMDAVecGetArray(da,b,&array);
1008: for (j=ys; j<ys+ym; j++) {
1009: for (i=xs; i<xs+xm; i++) {
1010: array[j][i] = PetscExpScalar(-((PetscReal)i*Hx)*((PetscReal)i*Hx)/user->nu)*PetscExpScalar(-((PetscReal)j*Hy)*((PetscReal)j*Hy)/user->nu)*Hx*Hy;
1011: }
1012: }
1013: DMDAVecRestoreArray(da, b, &array);
1014: VecAssemblyBegin(b);
1015: VecAssemblyEnd(b);
1017: /* force right hand side to be consistent for singular matrix */
1018: /* note this is really a hack, normally the model would provide you with a consistent right handside */
1019: if (user->bcType == NEUMANN) {
1020: MatNullSpace nullspace;
1022: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
1023: MatNullSpaceRemove(nullspace,b);
1024: MatNullSpaceDestroy(&nullspace);
1025: }
1026: return 0;
1027: }
1029: PetscErrorCode ComputeRho(PetscInt i, PetscInt j, PetscInt mx, PetscInt my, PetscReal centerRho, PetscReal *rho)
1030: {
1032: if ((i > mx/3.0) && (i < 2.0*mx/3.0) && (j > my/3.0) && (j < 2.0*my/3.0)) {
1033: *rho = centerRho;
1034: } else {
1035: *rho = 1.0;
1036: }
1037: return 0;
1038: }
1040: PetscErrorCode ComputeMatrix_DMDA(DM da,Mat J,Mat jac,void *ctx)
1041: {
1042: UserContext *user = (UserContext*)ctx;
1043: PetscReal centerRho;
1044: PetscInt i,j,mx,my,xm,ym,xs,ys;
1045: PetscScalar v[5];
1046: PetscReal Hx,Hy,HydHx,HxdHy,rho;
1047: MatStencil row, col[5];
1048: PetscBool isda = PETSC_FALSE;
1051: PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);
1053: MatZeroEntries(jac);
1054: centerRho = user->rho;
1055: DMDAGetInfo(da,NULL,&mx,&my,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1056: Hx = 1.0 / (PetscReal)(mx-1);
1057: Hy = 1.0 / (PetscReal)(my-1);
1058: HxdHy = Hx/Hy;
1059: HydHx = Hy/Hx;
1060: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
1061: for (j=ys; j<ys+ym; j++) {
1062: for (i=xs; i<xs+xm; i++) {
1063: row.i = i; row.j = j;
1064: ComputeRho(i, j, mx, my, centerRho, &rho);
1065: if (i==0 || j==0 || i==mx-1 || j==my-1) {
1066: if (user->bcType == DIRICHLET) {
1067: v[0] = 2.0*rho*(HxdHy + HydHx);
1068: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
1069: } else if (user->bcType == NEUMANN) {
1070: PetscInt numx = 0, numy = 0, num = 0;
1071: if (j!=0) {
1072: v[num] = -rho*HxdHy; col[num].i = i; col[num].j = j-1;
1073: numy++; num++;
1074: }
1075: if (i!=0) {
1076: v[num] = -rho*HydHx; col[num].i = i-1; col[num].j = j;
1077: numx++; num++;
1078: }
1079: if (i!=mx-1) {
1080: v[num] = -rho*HydHx; col[num].i = i+1; col[num].j = j;
1081: numx++; num++;
1082: }
1083: if (j!=my-1) {
1084: v[num] = -rho*HxdHy; col[num].i = i; col[num].j = j+1;
1085: numy++; num++;
1086: }
1087: v[num] = numx*rho*HydHx + numy*rho*HxdHy; col[num].i = i; col[num].j = j;
1088: num++;
1089: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
1090: }
1091: } else {
1092: v[0] = -rho*HxdHy; col[0].i = i; col[0].j = j-1;
1093: v[1] = -rho*HydHx; col[1].i = i-1; col[1].j = j;
1094: v[2] = 2.0*rho*(HxdHy + HydHx); col[2].i = i; col[2].j = j;
1095: v[3] = -rho*HydHx; col[3].i = i+1; col[3].j = j;
1096: v[4] = -rho*HxdHy; col[4].i = i; col[4].j = j+1;
1097: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
1098: }
1099: }
1100: }
1101: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
1102: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
1103: return 0;
1104: }
1106: PetscErrorCode ComputeMatrix_ShellDA(KSP ksp,Mat J,Mat jac,void *ctx)
1107: {
1108: DM dm,da;
1110: KSPGetDM(ksp,&dm);
1111: DMShellGetContext(dm,&da);
1112: ComputeMatrix_DMDA(da,J,jac,ctx);
1113: return 0;
1114: }
1116: /*TEST
1118: test:
1119: suffix: basic_dirichlet
1120: nsize: 4
1121: 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
1123: test:
1124: suffix: basic_neumann
1125: nsize: 4
1126: requires: !single
1127: 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
1129: test:
1130: suffix: mg_2lv_2mg
1131: nsize: 6
1132: 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
1134: test:
1135: suffix: mg_3lv_2mg
1136: nsize: 4
1137: 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
1139: test:
1140: suffix: mg_3lv_2mg_customcommsize
1141: nsize: 12
1142: 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
1144: TEST*/