Actual source code: wb.c

petsc-3.8.3 2017-12-09
Report Typos and Errors

  2:  #include <petscdmda.h>
  3:  #include <petsc/private/pcmgimpl.h>
  4:  #include <petscctable.h>

  6: typedef struct {
  7:   PCExoticType type;
  8:   Mat          P;            /* the constructed interpolation matrix */
  9:   PetscBool    directSolve;  /* use direct LU factorization to construct interpolation */
 10:   KSP          ksp;
 11: } PC_Exotic;

 13: const char *const PCExoticTypes[] = {"face","wirebasket","PCExoticType","PC_Exotic",0};


 16: /*
 17:       DMDAGetWireBasketInterpolation - Gets the interpolation for a wirebasket based coarse space

 19: */
 20: PetscErrorCode DMDAGetWireBasketInterpolation(DM da,PC_Exotic *exotic,Mat Aglobal,MatReuse reuse,Mat *P)
 21: {
 22:   PetscErrorCode         ierr;
 23:   PetscInt               dim,i,j,k,m,n,p,dof,Nint,Nface,Nwire,Nsurf,*Iint,*Isurf,cint = 0,csurf = 0,istart,jstart,kstart,*II,N,c = 0;
 24:   PetscInt               mwidth,nwidth,pwidth,cnt,mp,np,pp,Ntotal,gl[26],*globals,Ng,*IIint,*IIsurf,Nt;
 25:   Mat                    Xint, Xsurf,Xint_tmp;
 26:   IS                     isint,issurf,is,row,col;
 27:   ISLocalToGlobalMapping ltg;
 28:   MPI_Comm               comm;
 29:   Mat                    A,Aii,Ais,Asi,*Aholder,iAii;
 30:   MatFactorInfo          info;
 31:   PetscScalar            *xsurf,*xint;
 32: #if defined(PETSC_USE_DEBUG_foo)
 33:   PetscScalar            tmp;
 34: #endif
 35:   PetscTable             ht;

 38:   DMDAGetInfo(da,&dim,0,0,0,&mp,&np,&pp,&dof,0,0,0,0,0);
 39:   if (dof != 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only for single field problems");
 40:   if (dim != 3) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only coded for 3d problems");
 41:   DMDAGetCorners(da,0,0,0,&m,&n,&p);
 42:   DMDAGetGhostCorners(da,&istart,&jstart,&kstart,&mwidth,&nwidth,&pwidth);
 43:   istart = istart ? -1 : 0;
 44:   jstart = jstart ? -1 : 0;
 45:   kstart = kstart ? -1 : 0;

 47:   /*
 48:     the columns of P are the interpolation of each coarse grid point (one for each vertex and edge)
 49:     to all the local degrees of freedom (this includes the vertices, edges and faces).

 51:     Xint are the subset of the interpolation into the interior

 53:     Xface are the interpolation onto faces but not into the interior

 55:     Xsurf are the interpolation onto the vertices and edges (the surfbasket)
 56:                                       Xint
 57:     Symbolically one could write P = (Xface) after interchanging the rows to match the natural ordering on the domain
 58:                                       Xsurf
 59:   */
 60:   N     = (m - istart)*(n - jstart)*(p - kstart);
 61:   Nint  = (m-2-istart)*(n-2-jstart)*(p-2-kstart);
 62:   Nface = 2*((m-2-istart)*(n-2-jstart) + (m-2-istart)*(p-2-kstart) + (n-2-jstart)*(p-2-kstart));
 63:   Nwire = 4*((m-2-istart) + (n-2-jstart) + (p-2-kstart)) + 8;
 64:   Nsurf = Nface + Nwire;
 65:   MatCreateSeqDense(MPI_COMM_SELF,Nint,26,NULL,&Xint);
 66:   MatCreateSeqDense(MPI_COMM_SELF,Nsurf,26,NULL,&Xsurf);
 67:   MatDenseGetArray(Xsurf,&xsurf);

 69:   /*
 70:      Require that all 12 edges and 6 faces have at least one grid point. Otherwise some of the columns of
 71:      Xsurf will be all zero (thus making the coarse matrix singular).
 72:   */
 73:   if (m-istart < 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of grid points per process in X direction must be at least 3");
 74:   if (n-jstart < 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of grid points per process in Y direction must be at least 3");
 75:   if (p-kstart < 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of grid points per process in Z direction must be at least 3");

 77:   cnt = 0;

 79:   xsurf[cnt++] = 1;
 80:   for (i=1; i<m-istart-1; i++) xsurf[cnt++ + Nsurf] = 1;
 81:   xsurf[cnt++ + 2*Nsurf] = 1;

 83:   for (j=1; j<n-1-jstart; j++) {
 84:     xsurf[cnt++ + 3*Nsurf] = 1;
 85:     for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 4*Nsurf] = 1;
 86:     xsurf[cnt++ + 5*Nsurf] = 1;
 87:   }

 89:   xsurf[cnt++ + 6*Nsurf] = 1;
 90:   for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 7*Nsurf] = 1;
 91:   xsurf[cnt++ + 8*Nsurf] = 1;

 93:   for (k=1; k<p-1-kstart; k++) {
 94:     xsurf[cnt++ + 9*Nsurf] = 1;
 95:     for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 10*Nsurf] = 1;
 96:     xsurf[cnt++ + 11*Nsurf] = 1;

 98:     for (j=1; j<n-1-jstart; j++) {
 99:       xsurf[cnt++ + 12*Nsurf] = 1;
100:       /* these are the interior nodes */
101:       xsurf[cnt++ + 13*Nsurf] = 1;
102:     }

104:     xsurf[cnt++ + 14*Nsurf] = 1;
105:     for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 15*Nsurf] = 1;
106:     xsurf[cnt++ + 16*Nsurf] = 1;
107:   }

109:   xsurf[cnt++ + 17*Nsurf] = 1;
110:   for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 18*Nsurf] = 1;
111:   xsurf[cnt++ + 19*Nsurf] = 1;

113:   for (j=1;j<n-1-jstart;j++) {
114:     xsurf[cnt++ + 20*Nsurf] = 1;
115:     for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 21*Nsurf] = 1;
116:     xsurf[cnt++ + 22*Nsurf] = 1;
117:   }

119:   xsurf[cnt++ + 23*Nsurf] = 1;
120:   for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 24*Nsurf] = 1;
121:   xsurf[cnt++ + 25*Nsurf] = 1;


124:   /* interpolations only sum to 1 when using direct solver */
125: #if defined(PETSC_USE_DEBUG_foo)
126:   for (i=0; i<Nsurf; i++) {
127:     tmp = 0.0;
128:     for (j=0; j<26; j++) tmp += xsurf[i+j*Nsurf];
129:     if (PetscAbsScalar(tmp-1.0) > 1.e-10) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong Xsurf interpolation at i %D value %g",i,(double)PetscAbsScalar(tmp));
130:   }
131: #endif
132:   MatDenseRestoreArray(Xsurf,&xsurf);
133:   /* MatView(Xsurf,PETSC_VIEWER_STDOUT_WORLD);*/


136:   /*
137:        I are the indices for all the needed vertices (in global numbering)
138:        Iint are the indices for the interior values, I surf for the surface values
139:             (This is just for the part of the global matrix obtained with MatCreateSubMatrix(), it
140:              is NOT the local DMDA ordering.)
141:        IIint and IIsurf are the same as the Iint, Isurf except they are in the global numbering
142:   */
143: #define Endpoint(a,start,b) (a == 0 || a == (b-1-start))
144:   PetscMalloc3(N,&II,Nint,&Iint,Nsurf,&Isurf);
145:   PetscMalloc2(Nint,&IIint,Nsurf,&IIsurf);
146:   for (k=0; k<p-kstart; k++) {
147:     for (j=0; j<n-jstart; j++) {
148:       for (i=0; i<m-istart; i++) {
149:         II[c++] = i + j*mwidth + k*mwidth*nwidth;

151:         if (!Endpoint(i,istart,m) && !Endpoint(j,jstart,n) && !Endpoint(k,kstart,p)) {
152:           IIint[cint]  = i + j*mwidth + k*mwidth*nwidth;
153:           Iint[cint++] = i + j*(m-istart) + k*(m-istart)*(n-jstart);
154:         } else {
155:           IIsurf[csurf]  = i + j*mwidth + k*mwidth*nwidth;
156:           Isurf[csurf++] = i + j*(m-istart) + k*(m-istart)*(n-jstart);
157:         }
158:       }
159:     }
160:   }
161:   if (c != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"c != N");
162:   if (cint != Nint) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"cint != Nint");
163:   if (csurf != Nsurf) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"csurf != Nsurf");
164:   DMGetLocalToGlobalMapping(da,&ltg);
165:   ISLocalToGlobalMappingApply(ltg,N,II,II);
166:   ISLocalToGlobalMappingApply(ltg,Nint,IIint,IIint);
167:   ISLocalToGlobalMappingApply(ltg,Nsurf,IIsurf,IIsurf);
168:   PetscObjectGetComm((PetscObject)da,&comm);
169:   ISCreateGeneral(comm,N,II,PETSC_COPY_VALUES,&is);
170:   ISCreateGeneral(PETSC_COMM_SELF,Nint,Iint,PETSC_COPY_VALUES,&isint);
171:   ISCreateGeneral(PETSC_COMM_SELF,Nsurf,Isurf,PETSC_COPY_VALUES,&issurf);
172:   PetscFree3(II,Iint,Isurf);

174:   MatCreateSubMatrices(Aglobal,1,&is,&is,MAT_INITIAL_MATRIX,&Aholder);
175:   A    = *Aholder;
176:   PetscFree(Aholder);

178:   MatCreateSubMatrix(A,isint,isint,MAT_INITIAL_MATRIX,&Aii);
179:   MatCreateSubMatrix(A,isint,issurf,MAT_INITIAL_MATRIX,&Ais);
180:   MatCreateSubMatrix(A,issurf,isint,MAT_INITIAL_MATRIX,&Asi);

182:   /*
183:      Solve for the interpolation onto the interior Xint
184:   */
185:   MatMatMult(Ais,Xsurf,MAT_INITIAL_MATRIX,PETSC_DETERMINE,&Xint_tmp);
186:   MatScale(Xint_tmp,-1.0);
187:   if (exotic->directSolve) {
188:     MatGetFactor(Aii,MATSOLVERPETSC,MAT_FACTOR_LU,&iAii);
189:     MatFactorInfoInitialize(&info);
190:     MatGetOrdering(Aii,MATORDERINGND,&row,&col);
191:     MatLUFactorSymbolic(iAii,Aii,row,col,&info);
192:     ISDestroy(&row);
193:     ISDestroy(&col);
194:     MatLUFactorNumeric(iAii,Aii,&info);
195:     MatMatSolve(iAii,Xint_tmp,Xint);
196:     MatDestroy(&iAii);
197:   } else {
198:     Vec         b,x;
199:     PetscScalar *xint_tmp;

201:     MatDenseGetArray(Xint,&xint);
202:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,Nint,0,&x);
203:     MatDenseGetArray(Xint_tmp,&xint_tmp);
204:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,Nint,0,&b);
205:     KSPSetOperators(exotic->ksp,Aii,Aii);
206:     for (i=0; i<26; i++) {
207:       VecPlaceArray(x,xint+i*Nint);
208:       VecPlaceArray(b,xint_tmp+i*Nint);
209:       KSPSolve(exotic->ksp,b,x);
210:       VecResetArray(x);
211:       VecResetArray(b);
212:     }
213:     MatDenseRestoreArray(Xint,&xint);
214:     MatDenseRestoreArray(Xint_tmp,&xint_tmp);
215:     VecDestroy(&x);
216:     VecDestroy(&b);
217:   }
218:   MatDestroy(&Xint_tmp);

220: #if defined(PETSC_USE_DEBUG_foo)
221:   MatDenseGetArray(Xint,&xint);
222:   for (i=0; i<Nint; i++) {
223:     tmp = 0.0;
224:     for (j=0; j<26; j++) tmp += xint[i+j*Nint];

226:     if (PetscAbsScalar(tmp-1.0) > 1.e-10) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong Xint interpolation at i %D value %g",i,(double)PetscAbsScalar(tmp));
227:   }
228:   MatDenseRestoreArray(Xint,&xint);
229:   /* ierr =MatView(Xint,PETSC_VIEWER_STDOUT_WORLD); */
230: #endif


233:   /*         total vertices             total faces                                  total edges */
234:   Ntotal = (mp + 1)*(np + 1)*(pp + 1) + mp*np*(pp+1) + mp*pp*(np+1) + np*pp*(mp+1) + mp*(np+1)*(pp+1) + np*(mp+1)*(pp+1) +  pp*(mp+1)*(np+1);

236:   /*
237:       For each vertex, edge, face on process (in the same orderings as used above) determine its local number including ghost points
238:   */
239:   cnt = 0;

241:   gl[cnt++] = 0;  { gl[cnt++] = 1;} gl[cnt++] = m-istart-1;
242:   { gl[cnt++] = mwidth;  { gl[cnt++] = mwidth+1;} gl[cnt++] = mwidth + m-istart-1;}
243:   gl[cnt++] = mwidth*(n-jstart-1);  { gl[cnt++] = mwidth*(n-jstart-1)+1;} gl[cnt++] = mwidth*(n-jstart-1) + m-istart-1;
244:   {
245:     gl[cnt++] = mwidth*nwidth;  { gl[cnt++] = mwidth*nwidth+1;}  gl[cnt++] = mwidth*nwidth+ m-istart-1;
246:     { gl[cnt++] = mwidth*nwidth + mwidth; /* these are the interior nodes */ gl[cnt++] = mwidth*nwidth + mwidth+m-istart-1;}
247:     gl[cnt++] = mwidth*nwidth+ mwidth*(n-jstart-1);   { gl[cnt++] = mwidth*nwidth+mwidth*(n-jstart-1)+1;} gl[cnt++] = mwidth*nwidth+mwidth*(n-jstart-1) + m-istart-1;
248:   }
249:   gl[cnt++] = mwidth*nwidth*(p-kstart-1); { gl[cnt++] = mwidth*nwidth*(p-kstart-1)+1;} gl[cnt++] = mwidth*nwidth*(p-kstart-1) +  m-istart-1;
250:   { gl[cnt++] = mwidth*nwidth*(p-kstart-1) + mwidth;   { gl[cnt++] = mwidth*nwidth*(p-kstart-1) + mwidth+1;} gl[cnt++] = mwidth*nwidth*(p-kstart-1)+mwidth+m-istart-1;}
251:   gl[cnt++] = mwidth*nwidth*(p-kstart-1) +  mwidth*(n-jstart-1);  { gl[cnt++] = mwidth*nwidth*(p-kstart-1)+ mwidth*(n-jstart-1)+1;} gl[cnt++] = mwidth*nwidth*(p-kstart-1) + mwidth*(n-jstart-1) + m-istart-1;

253:   /* PetscIntView(26,gl,PETSC_VIEWER_STDOUT_WORLD); */
254:   /* convert that to global numbering and get them on all processes */
255:   ISLocalToGlobalMappingApply(ltg,26,gl,gl);
256:   /* PetscIntView(26,gl,PETSC_VIEWER_STDOUT_WORLD); */
257:   PetscMalloc1(26*mp*np*pp,&globals);
258:   MPI_Allgather(gl,26,MPIU_INT,globals,26,MPIU_INT,PetscObjectComm((PetscObject)da));

260:   /* Number the coarse grid points from 0 to Ntotal */
261:   MatGetSize(Aglobal,&Nt,NULL);
262:   PetscTableCreate(Ntotal/3,Nt+1,&ht);
263:   for (i=0; i<26*mp*np*pp; i++) {
264:     PetscTableAddCount(ht,globals[i]+1);
265:   }
266:   PetscTableGetCount(ht,&cnt);
267:   if (cnt != Ntotal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Hash table size %D not equal to total number coarse grid points %D",cnt,Ntotal);
268:   PetscFree(globals);
269:   for (i=0; i<26; i++) {
270:     PetscTableFind(ht,gl[i]+1,&gl[i]);
271:     gl[i]--;
272:   }
273:   PetscTableDestroy(&ht);
274:   /* PetscIntView(26,gl,PETSC_VIEWER_STDOUT_WORLD); */

276:   /* construct global interpolation matrix */
277:   MatGetLocalSize(Aglobal,&Ng,NULL);
278:   if (reuse == MAT_INITIAL_MATRIX) {
279:     MatCreateAIJ(PetscObjectComm((PetscObject)da),Ng,PETSC_DECIDE,PETSC_DECIDE,Ntotal,Nint+Nsurf,NULL,Nint+Nsurf,NULL,P);
280:   } else {
281:     MatZeroEntries(*P);
282:   }
283:   MatSetOption(*P,MAT_ROW_ORIENTED,PETSC_FALSE);
284:   MatDenseGetArray(Xint,&xint);
285:   MatSetValues(*P,Nint,IIint,26,gl,xint,INSERT_VALUES);
286:   MatDenseRestoreArray(Xint,&xint);
287:   MatDenseGetArray(Xsurf,&xsurf);
288:   MatSetValues(*P,Nsurf,IIsurf,26,gl,xsurf,INSERT_VALUES);
289:   MatDenseRestoreArray(Xsurf,&xsurf);
290:   MatAssemblyBegin(*P,MAT_FINAL_ASSEMBLY);
291:   MatAssemblyEnd(*P,MAT_FINAL_ASSEMBLY);
292:   PetscFree2(IIint,IIsurf);

294: #if defined(PETSC_USE_DEBUG_foo)
295:   {
296:     Vec         x,y;
297:     PetscScalar *yy;
298:     VecCreateMPI(PetscObjectComm((PetscObject)da),Ng,PETSC_DETERMINE,&y);
299:     VecCreateMPI(PetscObjectComm((PetscObject)da),PETSC_DETERMINE,Ntotal,&x);
300:     VecSet(x,1.0);
301:     MatMult(*P,x,y);
302:     VecGetArray(y,&yy);
303:     for (i=0; i<Ng; i++) {
304:       if (PetscAbsScalar(yy[i]-1.0) > 1.e-10) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong p interpolation at i %D value %g",i,(double)PetscAbsScalar(yy[i]));
305:     }
306:     VecRestoreArray(y,&yy);
307:     VecDestroy(x);
308:     VecDestroy(y);
309:   }
310: #endif

312:   MatDestroy(&Aii);
313:   MatDestroy(&Ais);
314:   MatDestroy(&Asi);
315:   MatDestroy(&A);
316:   ISDestroy(&is);
317:   ISDestroy(&isint);
318:   ISDestroy(&issurf);
319:   MatDestroy(&Xint);
320:   MatDestroy(&Xsurf);
321:   return(0);
322: }

324: /*
325:       DMDAGetFaceInterpolation - Gets the interpolation for a face based coarse space

327: */
328: PetscErrorCode DMDAGetFaceInterpolation(DM da,PC_Exotic *exotic,Mat Aglobal,MatReuse reuse,Mat *P)
329: {
330:   PetscErrorCode         ierr;
331:   PetscInt               dim,i,j,k,m,n,p,dof,Nint,Nface,Nwire,Nsurf,*Iint,*Isurf,cint = 0,csurf = 0,istart,jstart,kstart,*II,N,c = 0;
332:   PetscInt               mwidth,nwidth,pwidth,cnt,mp,np,pp,Ntotal,gl[6],*globals,Ng,*IIint,*IIsurf,Nt;
333:   Mat                    Xint, Xsurf,Xint_tmp;
334:   IS                     isint,issurf,is,row,col;
335:   ISLocalToGlobalMapping ltg;
336:   MPI_Comm               comm;
337:   Mat                    A,Aii,Ais,Asi,*Aholder,iAii;
338:   MatFactorInfo          info;
339:   PetscScalar            *xsurf,*xint;
340: #if defined(PETSC_USE_DEBUG_foo)
341:   PetscScalar            tmp;
342: #endif
343:   PetscTable             ht;

346:   DMDAGetInfo(da,&dim,0,0,0,&mp,&np,&pp,&dof,0,0,0,0,0);
347:   if (dof != 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only for single field problems");
348:   if (dim != 3) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only coded for 3d problems");
349:   DMDAGetCorners(da,0,0,0,&m,&n,&p);
350:   DMDAGetGhostCorners(da,&istart,&jstart,&kstart,&mwidth,&nwidth,&pwidth);
351:   istart = istart ? -1 : 0;
352:   jstart = jstart ? -1 : 0;
353:   kstart = kstart ? -1 : 0;

355:   /*
356:     the columns of P are the interpolation of each coarse grid point (one for each vertex and edge)
357:     to all the local degrees of freedom (this includes the vertices, edges and faces).

359:     Xint are the subset of the interpolation into the interior

361:     Xface are the interpolation onto faces but not into the interior

363:     Xsurf are the interpolation onto the vertices and edges (the surfbasket)
364:                                       Xint
365:     Symbolically one could write P = (Xface) after interchanging the rows to match the natural ordering on the domain
366:                                       Xsurf
367:   */
368:   N     = (m - istart)*(n - jstart)*(p - kstart);
369:   Nint  = (m-2-istart)*(n-2-jstart)*(p-2-kstart);
370:   Nface = 2*((m-2-istart)*(n-2-jstart) + (m-2-istart)*(p-2-kstart) + (n-2-jstart)*(p-2-kstart));
371:   Nwire = 4*((m-2-istart) + (n-2-jstart) + (p-2-kstart)) + 8;
372:   Nsurf = Nface + Nwire;
373:   MatCreateSeqDense(MPI_COMM_SELF,Nint,6,NULL,&Xint);
374:   MatCreateSeqDense(MPI_COMM_SELF,Nsurf,6,NULL,&Xsurf);
375:   MatDenseGetArray(Xsurf,&xsurf);

377:   /*
378:      Require that all 12 edges and 6 faces have at least one grid point. Otherwise some of the columns of
379:      Xsurf will be all zero (thus making the coarse matrix singular).
380:   */
381:   if (m-istart < 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of grid points per process in X direction must be at least 3");
382:   if (n-jstart < 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of grid points per process in Y direction must be at least 3");
383:   if (p-kstart < 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Number of grid points per process in Z direction must be at least 3");

385:   cnt = 0;
386:   for (j=1; j<n-1-jstart; j++) {
387:     for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 0*Nsurf] = 1;
388:   }

390:   for (k=1; k<p-1-kstart; k++) {
391:     for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 1*Nsurf] = 1;
392:     for (j=1; j<n-1-jstart; j++) {
393:       xsurf[cnt++ + 2*Nsurf] = 1;
394:       /* these are the interior nodes */
395:       xsurf[cnt++ + 3*Nsurf] = 1;
396:     }
397:     for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 4*Nsurf] = 1;
398:   }
399:   for (j=1;j<n-1-jstart;j++) {
400:     for (i=1; i<m-istart-1; i++) xsurf[cnt++ + 5*Nsurf] = 1;
401:   }

403: #if defined(PETSC_USE_DEBUG_foo)
404:   for (i=0; i<Nsurf; i++) {
405:     tmp = 0.0;
406:     for (j=0; j<6; j++) tmp += xsurf[i+j*Nsurf];

408:     if (PetscAbsScalar(tmp-1.0) > 1.e-10) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong Xsurf interpolation at i %D value %g",i,(double)PetscAbsScalar(tmp));
409:   }
410: #endif
411:   MatDenseRestoreArray(Xsurf,&xsurf);
412:   /* MatView(Xsurf,PETSC_VIEWER_STDOUT_WORLD);*/


415:   /*
416:        I are the indices for all the needed vertices (in global numbering)
417:        Iint are the indices for the interior values, I surf for the surface values
418:             (This is just for the part of the global matrix obtained with MatCreateSubMatrix(), it
419:              is NOT the local DMDA ordering.)
420:        IIint and IIsurf are the same as the Iint, Isurf except they are in the global numbering
421:   */
422: #define Endpoint(a,start,b) (a == 0 || a == (b-1-start))
423:   PetscMalloc3(N,&II,Nint,&Iint,Nsurf,&Isurf);
424:   PetscMalloc2(Nint,&IIint,Nsurf,&IIsurf);
425:   for (k=0; k<p-kstart; k++) {
426:     for (j=0; j<n-jstart; j++) {
427:       for (i=0; i<m-istart; i++) {
428:         II[c++] = i + j*mwidth + k*mwidth*nwidth;

430:         if (!Endpoint(i,istart,m) && !Endpoint(j,jstart,n) && !Endpoint(k,kstart,p)) {
431:           IIint[cint]  = i + j*mwidth + k*mwidth*nwidth;
432:           Iint[cint++] = i + j*(m-istart) + k*(m-istart)*(n-jstart);
433:         } else {
434:           IIsurf[csurf]  = i + j*mwidth + k*mwidth*nwidth;
435:           Isurf[csurf++] = i + j*(m-istart) + k*(m-istart)*(n-jstart);
436:         }
437:       }
438:     }
439:   }
440:   if (c != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"c != N");
441:   if (cint != Nint) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"cint != Nint");
442:   if (csurf != Nsurf) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"csurf != Nsurf");
443:   DMGetLocalToGlobalMapping(da,&ltg);
444:   ISLocalToGlobalMappingApply(ltg,N,II,II);
445:   ISLocalToGlobalMappingApply(ltg,Nint,IIint,IIint);
446:   ISLocalToGlobalMappingApply(ltg,Nsurf,IIsurf,IIsurf);
447:   PetscObjectGetComm((PetscObject)da,&comm);
448:   ISCreateGeneral(comm,N,II,PETSC_COPY_VALUES,&is);
449:   ISCreateGeneral(PETSC_COMM_SELF,Nint,Iint,PETSC_COPY_VALUES,&isint);
450:   ISCreateGeneral(PETSC_COMM_SELF,Nsurf,Isurf,PETSC_COPY_VALUES,&issurf);
451:   PetscFree3(II,Iint,Isurf);

453:   ISSort(is);
454:   MatCreateSubMatrices(Aglobal,1,&is,&is,MAT_INITIAL_MATRIX,&Aholder);
455:   A    = *Aholder;
456:   PetscFree(Aholder);

458:   MatCreateSubMatrix(A,isint,isint,MAT_INITIAL_MATRIX,&Aii);
459:   MatCreateSubMatrix(A,isint,issurf,MAT_INITIAL_MATRIX,&Ais);
460:   MatCreateSubMatrix(A,issurf,isint,MAT_INITIAL_MATRIX,&Asi);

462:   /*
463:      Solve for the interpolation onto the interior Xint
464:   */
465:   MatMatMult(Ais,Xsurf,MAT_INITIAL_MATRIX,PETSC_DETERMINE,&Xint_tmp);
466:   MatScale(Xint_tmp,-1.0);

468:   if (exotic->directSolve) {
469:     MatGetFactor(Aii,MATSOLVERPETSC,MAT_FACTOR_LU,&iAii);
470:     MatFactorInfoInitialize(&info);
471:     MatGetOrdering(Aii,MATORDERINGND,&row,&col);
472:     MatLUFactorSymbolic(iAii,Aii,row,col,&info);
473:     ISDestroy(&row);
474:     ISDestroy(&col);
475:     MatLUFactorNumeric(iAii,Aii,&info);
476:     MatMatSolve(iAii,Xint_tmp,Xint);
477:     MatDestroy(&iAii);
478:   } else {
479:     Vec         b,x;
480:     PetscScalar *xint_tmp;

482:     MatDenseGetArray(Xint,&xint);
483:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,Nint,0,&x);
484:     MatDenseGetArray(Xint_tmp,&xint_tmp);
485:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,Nint,0,&b);
486:     KSPSetOperators(exotic->ksp,Aii,Aii);
487:     for (i=0; i<6; i++) {
488:       VecPlaceArray(x,xint+i*Nint);
489:       VecPlaceArray(b,xint_tmp+i*Nint);
490:       KSPSolve(exotic->ksp,b,x);
491:       VecResetArray(x);
492:       VecResetArray(b);
493:     }
494:     MatDenseRestoreArray(Xint,&xint);
495:     MatDenseRestoreArray(Xint_tmp,&xint_tmp);
496:     VecDestroy(&x);
497:     VecDestroy(&b);
498:   }
499:   MatDestroy(&Xint_tmp);

501: #if defined(PETSC_USE_DEBUG_foo)
502:   MatDenseGetArray(Xint,&xint);
503:   for (i=0; i<Nint; i++) {
504:     tmp = 0.0;
505:     for (j=0; j<6; j++) tmp += xint[i+j*Nint];

507:     if (PetscAbsScalar(tmp-1.0) > 1.e-10) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong Xint interpolation at i %D value %g",i,(double)PetscAbsScalar(tmp));
508:   }
509:   MatDenseRestoreArray(Xint,&xint);
510:   /* ierr =MatView(Xint,PETSC_VIEWER_STDOUT_WORLD); */
511: #endif


514:   /*         total faces    */
515:   Ntotal =  mp*np*(pp+1) + mp*pp*(np+1) + np*pp*(mp+1);

517:   /*
518:       For each vertex, edge, face on process (in the same orderings as used above) determine its local number including ghost points
519:   */
520:   cnt = 0;
521:   { gl[cnt++] = mwidth+1;}
522:   {
523:     { gl[cnt++] = mwidth*nwidth+1;}
524:     { gl[cnt++] = mwidth*nwidth + mwidth; /* these are the interior nodes */ gl[cnt++] = mwidth*nwidth + mwidth+m-istart-1;}
525:     { gl[cnt++] = mwidth*nwidth+mwidth*(n-jstart-1)+1;}
526:   }
527:   { gl[cnt++] = mwidth*nwidth*(p-kstart-1) + mwidth+1;}

529:   /* PetscIntView(6,gl,PETSC_VIEWER_STDOUT_WORLD); */
530:   /* convert that to global numbering and get them on all processes */
531:   ISLocalToGlobalMappingApply(ltg,6,gl,gl);
532:   /* PetscIntView(6,gl,PETSC_VIEWER_STDOUT_WORLD); */
533:   PetscMalloc1(6*mp*np*pp,&globals);
534:   MPI_Allgather(gl,6,MPIU_INT,globals,6,MPIU_INT,PetscObjectComm((PetscObject)da));

536:   /* Number the coarse grid points from 0 to Ntotal */
537:   MatGetSize(Aglobal,&Nt,NULL);
538:   PetscTableCreate(Ntotal/3,Nt+1,&ht);
539:   for (i=0; i<6*mp*np*pp; i++) {
540:     PetscTableAddCount(ht,globals[i]+1);
541:   }
542:   PetscTableGetCount(ht,&cnt);
543:   if (cnt != Ntotal) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Hash table size %D not equal to total number coarse grid points %D",cnt,Ntotal);
544:   PetscFree(globals);
545:   for (i=0; i<6; i++) {
546:     PetscTableFind(ht,gl[i]+1,&gl[i]);
547:     gl[i]--;
548:   }
549:   PetscTableDestroy(&ht);
550:   /* PetscIntView(6,gl,PETSC_VIEWER_STDOUT_WORLD); */

552:   /* construct global interpolation matrix */
553:   MatGetLocalSize(Aglobal,&Ng,NULL);
554:   if (reuse == MAT_INITIAL_MATRIX) {
555:     MatCreateAIJ(PetscObjectComm((PetscObject)da),Ng,PETSC_DECIDE,PETSC_DECIDE,Ntotal,Nint+Nsurf,NULL,Nint,NULL,P);
556:   } else {
557:     MatZeroEntries(*P);
558:   }
559:   MatSetOption(*P,MAT_ROW_ORIENTED,PETSC_FALSE);
560:   MatDenseGetArray(Xint,&xint);
561:   MatSetValues(*P,Nint,IIint,6,gl,xint,INSERT_VALUES);
562:   MatDenseRestoreArray(Xint,&xint);
563:   MatDenseGetArray(Xsurf,&xsurf);
564:   MatSetValues(*P,Nsurf,IIsurf,6,gl,xsurf,INSERT_VALUES);
565:   MatDenseRestoreArray(Xsurf,&xsurf);
566:   MatAssemblyBegin(*P,MAT_FINAL_ASSEMBLY);
567:   MatAssemblyEnd(*P,MAT_FINAL_ASSEMBLY);
568:   PetscFree2(IIint,IIsurf);


571: #if defined(PETSC_USE_DEBUG_foo)
572:   {
573:     Vec         x,y;
574:     PetscScalar *yy;
575:     VecCreateMPI(PetscObjectComm((PetscObject)da),Ng,PETSC_DETERMINE,&y);
576:     VecCreateMPI(PetscObjectComm((PetscObject)da),PETSC_DETERMINE,Ntotal,&x);
577:     VecSet(x,1.0);
578:     MatMult(*P,x,y);
579:     VecGetArray(y,&yy);
580:     for (i=0; i<Ng; i++) {
581:       if (PetscAbsScalar(yy[i]-1.0) > 1.e-10) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Wrong p interpolation at i %D value %g",i,(double)PetscAbsScalar(yy[i]));
582:     }
583:     VecRestoreArray(y,&yy);
584:     VecDestroy(x);
585:     VecDestroy(y);
586:   }
587: #endif

589:   MatDestroy(&Aii);
590:   MatDestroy(&Ais);
591:   MatDestroy(&Asi);
592:   MatDestroy(&A);
593:   ISDestroy(&is);
594:   ISDestroy(&isint);
595:   ISDestroy(&issurf);
596:   MatDestroy(&Xint);
597:   MatDestroy(&Xsurf);
598:   return(0);
599: }


602: /*@
603:    PCExoticSetType - Sets the type of coarse grid interpolation to use

605:    Logically Collective on PC

607:    Input Parameters:
608: +  pc - the preconditioner context
609: -  type - either PC_EXOTIC_FACE or PC_EXOTIC_WIREBASKET (defaults to face)

611:    Notes: The face based interpolation has 1 degree of freedom per face and ignores the
612:      edge and vertex values completely in the coarse problem. For any seven point
613:      stencil the interpolation of a constant on all faces into the interior is that constant.

615:      The wirebasket interpolation has 1 degree of freedom per vertex, per edge and
616:      per face. A constant on the subdomain boundary is interpolated as that constant
617:      in the interior of the domain.

619:      The coarse grid matrix is obtained via the Galerkin computation A_c = R A R^T, hence
620:      if A is nonsingular A_c is also nonsingular.

622:      Both interpolations are suitable for only scalar problems.

624:    Level: intermediate


627: .seealso: PCEXOTIC, PCExoticType()
628: @*/
629: PetscErrorCode  PCExoticSetType(PC pc,PCExoticType type)
630: {

636:   PetscTryMethod(pc,"PCExoticSetType_C",(PC,PCExoticType),(pc,type));
637:   return(0);
638: }

640: static PetscErrorCode  PCExoticSetType_Exotic(PC pc,PCExoticType type)
641: {
642:   PC_MG     *mg  = (PC_MG*)pc->data;
643:   PC_Exotic *ctx = (PC_Exotic*) mg->innerctx;

646:   ctx->type = type;
647:   return(0);
648: }

650: PetscErrorCode PCSetUp_Exotic(PC pc)
651: {
653:   Mat            A;
654:   PC_MG          *mg   = (PC_MG*)pc->data;
655:   PC_Exotic      *ex   = (PC_Exotic*) mg->innerctx;
656:   MatReuse       reuse = (ex->P) ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;

659:   if (!pc->dm) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Need to call PCSetDM() before using this PC");
660:   PCGetOperators(pc,NULL,&A);
661:   if (ex->type == PC_EXOTIC_FACE) {
662:     DMDAGetFaceInterpolation(pc->dm,ex,A,reuse,&ex->P);
663:   } else if (ex->type == PC_EXOTIC_WIREBASKET) {
664:     DMDAGetWireBasketInterpolation(pc->dm,ex,A,reuse,&ex->P);
665:   } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Unknown exotic coarse space %d",ex->type);
666:   PCMGSetInterpolation(pc,1,ex->P);
667:   /* if PC has attached DM we must remove it or the PCMG will use it to compute incorrect sized vectors and interpolations */
668:   PCSetDM(pc,NULL);
669:   PCSetUp_MG(pc);
670:   return(0);
671: }

673: PetscErrorCode PCDestroy_Exotic(PC pc)
674: {
676:   PC_MG          *mg  = (PC_MG*)pc->data;
677:   PC_Exotic      *ctx = (PC_Exotic*) mg->innerctx;

680:   MatDestroy(&ctx->P);
681:   KSPDestroy(&ctx->ksp);
682:   PetscFree(ctx);
683:   PCDestroy_MG(pc);
684:   return(0);
685: }

687: PetscErrorCode PCView_Exotic(PC pc,PetscViewer viewer)
688: {
689:   PC_MG          *mg = (PC_MG*)pc->data;
691:   PetscBool      iascii;
692:   PC_Exotic      *ctx = (PC_Exotic*) mg->innerctx;

695:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
696:   if (iascii) {
697:     PetscViewerASCIIPrintf(viewer,"    Exotic type = %s\n",PCExoticTypes[ctx->type]);
698:     if (ctx->directSolve) {
699:       PetscViewerASCIIPrintf(viewer,"      Using direct solver to construct interpolation\n");
700:     } else {
701:       PetscViewer sviewer;
702:       PetscMPIInt rank;

704:       PetscViewerASCIIPrintf(viewer,"      Using iterative solver to construct interpolation\n");
705:       PetscViewerASCIIPushTab(viewer);
706:       PetscViewerASCIIPushTab(viewer);  /* should not need to push this twice? */
707:       PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
708:       MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
709:       if (!rank) {
710:         KSPView(ctx->ksp,sviewer);
711:       }
712:       PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
713:       PetscViewerASCIIPopTab(viewer);
714:       PetscViewerASCIIPopTab(viewer);
715:     }
716:   }
717:   PCView_MG(pc,viewer);
718:   return(0);
719: }

721: PetscErrorCode PCSetFromOptions_Exotic(PetscOptionItems *PetscOptionsObject,PC pc)
722: {
724:   PetscBool      flg;
725:   PC_MG          *mg = (PC_MG*)pc->data;
726:   PCExoticType   mgctype;
727:   PC_Exotic      *ctx = (PC_Exotic*) mg->innerctx;

730:   PetscOptionsHead(PetscOptionsObject,"Exotic coarse space options");
731:   PetscOptionsEnum("-pc_exotic_type","face or wirebasket","PCExoticSetType",PCExoticTypes,(PetscEnum)ctx->type,(PetscEnum*)&mgctype,&flg);
732:   if (flg) {
733:     PCExoticSetType(pc,mgctype);
734:   }
735:   PetscOptionsBool("-pc_exotic_direct_solver","use direct solver to construct interpolation","None",ctx->directSolve,&ctx->directSolve,NULL);
736:   if (!ctx->directSolve) {
737:     if (!ctx->ksp) {
738:       const char *prefix;
739:       KSPCreate(PETSC_COMM_SELF,&ctx->ksp);
740:       KSPSetErrorIfNotConverged(ctx->ksp,pc->erroriffailure);
741:       PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)pc,1);
742:       PetscLogObjectParent((PetscObject)pc,(PetscObject)ctx->ksp);
743:       PCGetOptionsPrefix(pc,&prefix);
744:       KSPSetOptionsPrefix(ctx->ksp,prefix);
745:       KSPAppendOptionsPrefix(ctx->ksp,"exotic_");
746:     }
747:     KSPSetFromOptions(ctx->ksp);
748:   }
749:   PetscOptionsTail();
750:   return(0);
751: }


754: /*MC
755:      PCEXOTIC - Two level overlapping Schwarz preconditioner with exotic (non-standard) coarse grid spaces

757:      This uses the PCMG infrastructure restricted to two levels and the face and wirebasket based coarse
758:    grid spaces.

760:    Notes: By default this uses GMRES on the fine grid smoother so this should be used with KSPFGMRES or the smoother changed to not use GMRES

762:    References:
763: +  1. - These coarse grid spaces originate in the work of Bramble, Pasciak  and Schatz, "The Construction
764:    of Preconditioners for Elliptic Problems by Substructing IV", Mathematics of Computation, volume 53, 1989.
765: .  2. - They were generalized slightly in "Domain Decomposition Method for Linear Elasticity", Ph. D. thesis, Barry Smith,
766:    New York University, 1990. 
767: .  3. - They were then explored in great detail in Dryja, Smith, Widlund, "Schwarz Analysis
768:    of Iterative Substructuring Methods for Elliptic Problems in Three Dimensions, SIAM Journal on Numerical
769:    Analysis, volume 31. 1994. These were developed in the context of iterative substructuring preconditioners.
770: .  4. - They were then ingeniously applied as coarse grid spaces for overlapping Schwarz methods by Dohrmann and Widlund.
771:    They refer to them as GDSW (generalized Dryja, Smith, Widlund preconditioners). See, for example,
772:    Clark R. Dohrmann, Axel Klawonn, and Olof B. Widlund. Extending theory for domain decomposition algorithms to irregular subdomains. In Ulrich Langer, Marco
773:    Discacciati, David Keyes, Olof Widlund, and Walter Zulehner, editors, Proceedings
774:    of the 17th International Conference on Domain Decomposition Methods in
775:    Science and Engineering, held in Strobl, Austria, 2006, number 60 in
776:    Springer Verlag, Lecture Notes in Computational Science and Engineering, 2007.
777: .  5. -  Clark R. Dohrmann, Axel Klawonn, and Olof B. Widlund. A family of energy minimizing coarse spaces for overlapping Schwarz preconditioners. In Ulrich Langer,
778:    Marco Discacciati, David Keyes, Olof Widlund, and Walter Zulehner, editors, Proceedings
779:    of the 17th International Conference on Domain Decomposition Methods
780:    in Science and Engineering, held in Strobl, Austria, 2006, number 60 in
781:    Springer Verlag, Lecture Notes in Computational Science and Engineering, 2007
782: .  6. - Clark R. Dohrmann, Axel Klawonn, and Olof B. Widlund. Domain decomposition
783:    for less regular subdomains: Overlapping Schwarz in two dimensions. SIAM J.
784:    Numer. Anal., 46(4), 2008.
785: -  7. - Clark R. Dohrmann and Olof B. Widlund. An overlapping Schwarz
786:    algorithm for almost incompressible elasticity. Technical Report
787:    TR2008 912, Department of Computer Science, Courant Institute
788:    of Mathematical Sciences, New York University, May 2008. URL:

790:    Options Database: The usual PCMG options are supported, such as -mg_levels_pc_type <type> -mg_coarse_pc_type <type>
791:       -pc_mg_type <type>

793:    Level: advanced

795: .seealso:  PCMG, PCSetDM(), PCExoticType, PCExoticSetType()
796: M*/

798: PETSC_EXTERN PetscErrorCode PCCreate_Exotic(PC pc)
799: {
801:   PC_Exotic      *ex;
802:   PC_MG          *mg;

805:   /* if type was previously mg; must manually destroy it because call to PCSetType(pc,PCMG) will not destroy it */
806:   if (pc->ops->destroy) {
807:      (*pc->ops->destroy)(pc);
808:     pc->data = 0;
809:   }
810:   PetscFree(((PetscObject)pc)->type_name);
811:   ((PetscObject)pc)->type_name = 0;

813:   PCSetType(pc,PCMG);
814:   PCMGSetLevels(pc,2,NULL);
815:   PCMGSetGalerkin(pc,PC_MG_GALERKIN_PMAT);
816:   PetscNew(&ex); \
817:   ex->type     = PC_EXOTIC_FACE;
818:   mg           = (PC_MG*) pc->data;
819:   mg->innerctx = ex;


822:   pc->ops->setfromoptions = PCSetFromOptions_Exotic;
823:   pc->ops->view           = PCView_Exotic;
824:   pc->ops->destroy        = PCDestroy_Exotic;
825:   pc->ops->setup          = PCSetUp_Exotic;

827:   PetscObjectComposeFunction((PetscObject)pc,"PCExoticSetType_C",PCExoticSetType_Exotic);
828:   return(0);
829: }