Actual source code: ex43.c

petsc-3.3-p7 2013-05-11
  1: static char help[] =
  2:   "Solves the incompressible, variable viscosity stokes equation in 2d on the unit domain \n\
  3: using Q1Q1 elements, stabilized with Bochev's polynomial projection method. \n\
  4: The models defined utilise free slip boundary conditions on all sides. \n\
  5: Options: \n\
  6:      -mx : number elements in x-direciton \n\
  7:      -my : number elements in y-direciton \n\
  8:      -c_str : indicates the structure of the coefficients to use. \n\
  9:           -c_str 0 => Setup for an analytic solution with a vertical jump in viscosity. This problem is driven by the \n\
 10:                          forcing function f = ( 0, sin(n_z pi y )cos( pi x ). \n\
 11:                          Parameters: \n\
 12:                               -solcx_eta0 : the viscosity to the left of the interface \n\
 13:                               -solcx_eta1 : the viscosity to the right of the interface \n\
 14:                               -solcx_xc : the location of the interface \n\
 15:                               -solcx_nz : the wavenumber in the y direction \n\
 16:           -c_str 1 => Setup for a rectangular blob located in the origin of the domain. \n\
 17:                          Parameters: \n\
 18:                               -sinker_eta0 : the viscosity of the background fluid \n\
 19:                               -sinker_eta1 : the viscosity of the blob \n\
 20:                               -sinker_dx : the width of the blob \n\
 21:                               -sinker_dy : the width of the blob \n\
 22:           -c_str 2 => Setup for a circular blob located in the origin of the domain. \n\
 23:                          Parameters: \n\
 24:                               -sinker_eta0 : the viscosity of the background fluid \n\
 25:                               -sinker_eta1 : the viscosity of the blob \n\
 26:                               -sinker_r : radius of the blob \n\
 27:           -c_str 3 => Circular and rectangular inclusion\n\
 28:                          Parameters as for cases 1 and 2 (above)\n\
 29:      -use_gp_coords : evaluate the viscosity and the body force at the global coordinates of the quadrature points.\n\
 30:      By default, eta and the body force are evaulated at the element center and applied as a constant over the entire element.\n";

 32: /* Contributed by Dave May */

 34: #include <petscksp.h>
 35: #include <petscdmda.h>

 37: /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
 38: #include "ex43-solCx.h"

 40: static PetscErrorCode DMDABCApplyFreeSlip(DM,Mat,Vec);


 43: #define NSD            2 /* number of spatial dimensions */
 44: #define NODES_PER_EL   4 /* nodes per element */
 45: #define U_DOFS         2 /* degrees of freedom per velocity node */
 46: #define P_DOFS         1 /* degrees of freedom per pressure node */
 47: #define GAUSS_POINTS   4

 49: /* cell based evaluation */
 50: typedef struct {
 51:   PetscScalar eta,fx,fy;
 52: } Coefficients;

 54: /* Gauss point based evaluation 8+4+4+4 = 20 */
 55: typedef struct {
 56:   PetscScalar gp_coords[2*GAUSS_POINTS];
 57:   PetscScalar eta[GAUSS_POINTS];
 58:   PetscScalar fx[GAUSS_POINTS];
 59:   PetscScalar fy[GAUSS_POINTS];
 60: } GaussPointCoefficients;

 62: typedef struct {
 63:   PetscScalar u_dof;
 64:   PetscScalar v_dof;
 65:   PetscScalar p_dof;
 66: } StokesDOF;


 69: /*

 71: D = [ 2.eta   0   0   ]
 72: [   0   2.eta 0   ]
 73: [   0     0   eta ]

 75: B = [ d_dx   0   ]
 76: [  0    d_dy ]
 77: [ d_dy  d_dx ]

 79: */

 81: /* FEM routines */
 82: /*
 83: Element: Local basis function ordering
 84: 1-----2
 85: |     |
 86: |     |
 87: 0-----3
 88: */
 89: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
 90: {
 91:   PetscScalar xi  = _xi[0];
 92:   PetscScalar eta = _xi[1];

 94:   Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
 95:   Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
 96:   Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
 97:   Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
 98: }

100: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
101: {
102:   PetscScalar xi  = _xi[0];
103:   PetscScalar eta = _xi[1];

105:   GNi[0][0] = -0.25*(1.0-eta);
106:   GNi[0][1] = -0.25*(1.0+eta);
107:   GNi[0][2] =   0.25*(1.0+eta);
108:   GNi[0][3] =   0.25*(1.0-eta);

110:   GNi[1][0] = -0.25*(1.0-xi);
111:   GNi[1][1] =   0.25*(1.0-xi);
112:   GNi[1][2] =   0.25*(1.0+xi);
113:   GNi[1][3] = -0.25*(1.0+xi);
114: }

116: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
117: {
118:   PetscScalar J00,J01,J10,J11,J;
119:   PetscScalar iJ00,iJ01,iJ10,iJ11;
120:   PetscInt    i;

122:   J00 = J01 = J10 = J11 = 0.0;
123:   for (i = 0; i < NODES_PER_EL; i++) {
124:     PetscScalar cx = coords[ 2*i+0 ];
125:     PetscScalar cy = coords[ 2*i+1 ];

127:     J00 = J00+GNi[0][i]*cx;      /* J_xx = dx/dxi */
128:     J01 = J01+GNi[0][i]*cy;      /* J_xy = dy/dxi */
129:     J10 = J10+GNi[1][i]*cx;      /* J_yx = dx/deta */
130:     J11 = J11+GNi[1][i]*cy;      /* J_yy = dy/deta */
131:   }
132:   J = (J00*J11)-(J01*J10);

134:   iJ00 =  J11/J;
135:   iJ01 = -J01/J;
136:   iJ10 = -J10/J;
137:   iJ11 =  J00/J;


140:   for (i = 0; i < NODES_PER_EL; i++) {
141:     GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
142:     GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
143:   }

145:   if (det_J != NULL) {
146:     *det_J = J;
147:   }
148: }

150: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
151: {
152:   *ngp         = 4;
153:   gp_xi[0][0]  = -0.57735026919;gp_xi[0][1] = -0.57735026919;
154:   gp_xi[1][0]  = -0.57735026919;gp_xi[1][1] =  0.57735026919;
155:   gp_xi[2][0]  =  0.57735026919;gp_xi[2][1] =  0.57735026919;
156:   gp_xi[3][0]  =  0.57735026919;gp_xi[3][1] = -0.57735026919;
157:   gp_weight[0] = 1.0;
158:   gp_weight[1] = 1.0;
159:   gp_weight[2] = 1.0;
160:   gp_weight[3] = 1.0;
161: }


164: /* procs to the left claim the ghost node as their element */
167: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
168: {
169:   PetscInt m,n,p,M,N,P;
170:   PetscInt sx,sy,sz;

173:   DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
174:   DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);

176:   if (mxl != PETSC_NULL) {
177:     *mxl = m;
178:     if ((sx+m) == M) {  /* last proc */
179:       *mxl = m-1;
180:     }
181:   }
182:   if (myl != PETSC_NULL) {
183:     *myl = n;
184:     if ((sy+n) == N) {  /* last proc */
185:       *myl = n-1;
186:     }
187:   }
188:   if (mzl != PETSC_NULL) {
189:     *mzl = p;
190:     if ((sz+p) == P) {  /* last proc */
191:       *mzl = p-1;
192:     }
193:   }
194:   return(0);
195: }

199: static PetscErrorCode DMDAGetElementCorners(DM da,
200:                                           PetscInt *sx,PetscInt *sy,PetscInt *sz,
201:                                           PetscInt *mx,PetscInt *my,PetscInt *mz)
202: {
203:   PetscInt si,sj,sk;

206:   DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);

208:   *sx = si;
209:   if (si != 0) {
210:     *sx = si+1;
211:   }

213:   *sy = sj;
214:   if (sj != 0) {
215:     *sy = sj+1;
216:   }

218:   if (sk != PETSC_NULL) {
219:     *sz = sk;
220:     if (sk != 0) {
221:       *sz = sk+1;
222:     }
223:   }

225:   DMDAGetLocalElementSize(da,mx,my,mz);

227:   return(0);
228: }

230: /*
231: i,j are the element indices
232: The unknown is a vector quantity.
233: The s[].c is used to indicate the degree of freedom.
234: */
237: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
238: {
240:   /* velocity */
241:   /* node 0 */
242:   s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0;                         /* Vx0 */
243:   s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1;                         /* Vy0 */

245:   /* node 1 */
246:   s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0;                         /* Vx1 */
247:   s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1;                         /* Vy1 */

249:   /* node 2 */
250:   s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0;                         /* Vx2 */
251:   s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1;                         /* Vy2 */

253:   /* node 3 */
254:   s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0;                         /* Vx3 */
255:   s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1;                         /* Vy3 */


258:   /* pressure */
259:   s_p[0].i = i;s_p[0].j = j;s_p[0].c = 2;                         /* P0 */
260:   s_p[1].i = i;s_p[1].j = j+1;s_p[1].c = 2;                         /* P0 */
261:   s_p[2].i = i+1;s_p[2].j = j+1;s_p[2].c = 2;                         /* P1 */
262:   s_p[3].i = i+1;s_p[3].j = j;s_p[3].c = 2;                         /* P1 */
263:   return(0);
264: }

268: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
269: {
271:   PetscMPIInt    rank;
272:   PetscInt       proc_I,proc_J;
273:   PetscInt       cpu_x,cpu_y;
274:   PetscInt       local_mx,local_my;
275:   Vec            vlx,vly;
276:   PetscInt       *LX,*LY,i;
277:   PetscScalar    *_a;
278:   Vec            V_SEQ;
279:   VecScatter     ctx;

282:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

284:   DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);

286:   proc_J = rank/cpu_x;
287:   proc_I = rank-cpu_x*proc_J;

289:   PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);
290:   PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);

292:   DMDAGetLocalElementSize(da,&local_mx,&local_my,PETSC_NULL);
293:   VecCreate(PETSC_COMM_WORLD,&vlx);
294:   VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
295:   VecSetFromOptions(vlx);

297:   VecCreate(PETSC_COMM_WORLD,&vly);
298:   VecSetSizes(vly,PETSC_DECIDE,cpu_y);
299:   VecSetFromOptions(vly);

301:   VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
302:   VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
303:   VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
304:   VecAssemblyBegin(vly);VecAssemblyEnd(vly);



308:   VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
309:   VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
310:   VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
311:   VecGetArray(V_SEQ,&_a);
312:   for (i = 0; i < cpu_x; i++) {
313:     LX[i] = (PetscInt)PetscRealPart(_a[i]);
314:   }
315:   VecRestoreArray(V_SEQ,&_a);
316:   VecScatterDestroy(&ctx);
317:   VecDestroy(&V_SEQ);

319:   VecScatterCreateToAll(vly,&ctx,&V_SEQ);
320:   VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
321:   VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
322:   VecGetArray(V_SEQ,&_a);
323:   for (i = 0; i < cpu_y; i++) {
324:     LY[i] = (PetscInt)PetscRealPart(_a[i]);
325:   }
326:   VecRestoreArray(V_SEQ,&_a);
327:   VecScatterDestroy(&ctx);
328:   VecDestroy(&V_SEQ);



332:   *_lx = LX;
333:   *_ly = LY;

335:   VecDestroy(&vlx);
336:   VecDestroy(&vly);

338:   return(0);
339: }

343: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
344: {
345:   DM             cda;
346:   Vec            coords;
347:   DMDACoor2d       **_coords;
348:   PetscInt       si,sj,nx,ny,i,j;
349:   FILE           *fp;
350:   char           fname[PETSC_MAX_PATH_LEN];
351:   PetscMPIInt    rank;

355:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
356:   PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
357:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
358:   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");

360:   PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);

362:   DMDAGetCoordinateDA(da,&cda);
363:   DMDAGetGhostedCoordinates(da,&coords);
364:   DMDAVecGetArray(cda,coords,&_coords);
365:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
366:   for (j = sj; j < sj+ny-1; j++) {
367:     for (i = si; i < si+nx-1; i++) {
368:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
369:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i].x),PetscRealPart(_coords[j+1][i].y));
370:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i+1].x),PetscRealPart(_coords[j+1][i+1].y));
371:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i+1].x),PetscRealPart(_coords[j][i+1].y));
372:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
373:     }
374:   }
375:   DMDAVecRestoreArray(cda,coords,&_coords);

377:   PetscFClose(PETSC_COMM_SELF,fp);
378:   return(0);
379: }

383: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
384: {
385:   DM             cda;
386:   Vec            coords,local_fields;
387:   DMDACoor2d       **_coords;
388:   FILE           *fp;
389:   char           fname[PETSC_MAX_PATH_LEN];
390:   PetscMPIInt    rank;
391:   PetscInt       si,sj,nx,ny,i,j;
392:   PetscInt       n_dofs,d;
393:   PetscScalar    *_fields;

397:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
398:   PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
399:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
400:   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");

402:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
403:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
404:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
405:   for (d = 0; d < n_dofs; d++) {
406:     const char *field_name;
407:     DMDAGetFieldName(da,d,&field_name);
408:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
409:   }
410:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");


413:   DMDAGetCoordinateDA(da,&cda);
414:   DMDAGetGhostedCoordinates(da,&coords);
415:   DMDAVecGetArray(cda,coords,&_coords);
416:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

418:   DMCreateLocalVector(da,&local_fields);
419:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
420:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
421:   VecGetArray(local_fields,&_fields);


424:   for (j = sj; j < sj+ny; j++) {
425:     for (i = si; i < si+nx; i++) {
426:       PetscScalar coord_x,coord_y;
427:       PetscScalar field_d;

429:       coord_x = _coords[j][i].x;
430:       coord_y = _coords[j][i].y;

432:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
433:       for (d = 0; d < n_dofs; d++) {
434:         field_d = _fields[ n_dofs*((i-si)+(j-sj)*(nx))+d ];
435:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",PetscRealPart(field_d));
436:       }
437:       PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
438:     }
439:   }
440:   VecRestoreArray(local_fields,&_fields);
441:   VecDestroy(&local_fields);

443:   DMDAVecRestoreArray(cda,coords,&_coords);

445:   PetscFClose(PETSC_COMM_SELF,fp);
446:   return(0);
447: }

451: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
452: {
453:   DM                     cda;
454:   Vec                    local_fields;
455:   FILE                   *fp;
456:   char                   fname[PETSC_MAX_PATH_LEN];
457:   PetscMPIInt            rank;
458:   PetscInt               si,sj,nx,ny,i,j,p;
459:   PetscInt               n_dofs,d;
460:   GaussPointCoefficients **_coefficients;
461:   PetscErrorCode         ierr;

464:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
465:   PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
466:   PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
467:   if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");

469:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
470:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
471:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
472:   for (d = 0; d < n_dofs; d++) {
473:     const char *field_name;
474:     DMDAGetFieldName(da,d,&field_name);
475:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
476:   }
477:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");


480:   DMDAGetCoordinateDA(da,&cda);
481:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

483:   DMCreateLocalVector(da,&local_fields);
484:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
485:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
486:   DMDAVecGetArray(da,local_fields,&_coefficients);


489:   for (j = sj; j < sj+ny; j++) {
490:     for (i = si; i < si+nx; i++) {
491:       PetscScalar coord_x,coord_y;

493:       for (p = 0; p < GAUSS_POINTS; p++) {
494:         coord_x = _coefficients[j][i].gp_coords[2*p];
495:         coord_y = _coefficients[j][i].gp_coords[2*p+1];

497:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));

499:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e",PetscRealPart(_coefficients[j][i].eta[p]),PetscRealPart(_coefficients[j][i].fx[p]),PetscRealPart(_coefficients[j][i].fy[p]));
500:         PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
501:       }
502:     }
503:   }
504:   DMDAVecRestoreArray(da,local_fields,&_coefficients);
505:   VecDestroy(&local_fields);

507:   PetscFClose(PETSC_COMM_SELF,fp);
508:   return(0);
509: }


512: static PetscInt ASS_MAP_wIwDI_uJuDJ(
513:   PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,
514:   PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
515: {
516:   PetscInt ij;
517:   PetscInt r,c,nc;

519:   nc = u_NPE*u_dof;

521:   r = w_dof*wi+wd;
522:   c = u_dof*ui+ud;

524:   ij = r*nc+c;

526:   return ij;
527: }

529: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
530: {
531:   PetscInt    ngp;
532:   PetscScalar gp_xi[GAUSS_POINTS][2];
533:   PetscScalar gp_weight[GAUSS_POINTS];
534:   PetscInt    p,i,j,k;
535:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
536:   PetscScalar J_p,tildeD[3];
537:   PetscScalar B[3][U_DOFS*NODES_PER_EL];


540:   /* define quadrature rule */
541:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

543:   /* evaluate integral */
544:   for (p = 0; p < ngp; p++) {
545:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
546:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);

548:     for (i = 0; i < NODES_PER_EL; i++) {
549:       PetscScalar d_dx_i = GNx_p[0][i];
550:       PetscScalar d_dy_i = GNx_p[1][i];

552:       B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
553:       B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
554:       B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
555:     }


558:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
559:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
560:     tildeD[2] =       gp_weight[p]*J_p*eta[p];

562:     /* form Bt tildeD B */
563:     /*
564:     Ke_ij = Bt_ik . D_kl . B_lj
565:     = B_ki . D_kl . B_lj
566:     = B_ki . D_kk . B_kj
567:     */
568:     for (i = 0; i < 8; i++) {
569:       for (j = 0; j < 8; j++) {
570:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
571:           Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
572:         }
573:       }
574:     }
575:   }
576: }

578: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
579: {
580:   PetscInt    ngp;
581:   PetscScalar gp_xi[GAUSS_POINTS][2];
582:   PetscScalar gp_weight[GAUSS_POINTS];
583:   PetscInt    p,i,j,di;
584:   PetscScalar Ni_p[NODES_PER_EL];
585:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
586:   PetscScalar J_p,fac;


589:   /* define quadrature rule */
590:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

592:   /* evaluate integral */
593:   for (p = 0; p < ngp; p++) {
594:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
595:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
596:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
597:     fac = gp_weight[p]*J_p;

599:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
600:       for (di = 0; di < NSD; di++) { /* u dofs */
601:         for (j = 0; j < 4; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
602:           PetscInt IJ;
603:           /*     Ke[4*u_idx+j] = Ke[4*u_idx+j] - GNx_p[di][i] * Ni_p[j] * fac; */
604:           IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);

606:           Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
607:         }
608:       }
609:     }
610:   }
611: }

613: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
614: {
615:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
616:   PetscInt    i,j;
617:   PetscInt    nr_g,nc_g;

619:   PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
620:   FormGradientOperatorQ1(Ge,coords);

622:   nr_g = U_DOFS*NODES_PER_EL;
623:   nc_g = P_DOFS*NODES_PER_EL;

625:   for (i = 0; i < nr_g; i++) {
626:     for (j = 0; j < nc_g; j++) {
627:       De[nr_g*j+i] = Ge[nc_g*i+j];
628:     }
629:   }
630: }

632: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
633: {
634:   PetscInt    ngp;
635:   PetscScalar gp_xi[GAUSS_POINTS][2];
636:   PetscScalar gp_weight[GAUSS_POINTS];
637:   PetscInt    p,i,j;
638:   PetscScalar Ni_p[NODES_PER_EL];
639:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
640:   PetscScalar J_p,fac,eta_avg;


643:   /* define quadrature rule */
644:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

646:   /* evaluate integral */
647:   for (p = 0; p < ngp; p++) {
648:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
649:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
650:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
651:     fac = gp_weight[p]*J_p;

653:     for (i = 0; i < NODES_PER_EL; i++) {
654:       for (j = 0; j < NODES_PER_EL; j++) {
655:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
656:       }
657:     }
658:   }

660:   /* scale */
661:   eta_avg = 0.0;
662:   for (p = 0; p < ngp; p++) {
663:     eta_avg += eta[p];
664:   }
665:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
666:   fac     = 1.0/eta_avg;
667:   for (i = 0; i < NODES_PER_EL; i++) {
668:     for (j = 0; j < NODES_PER_EL; j++) {
669:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
670:     }
671:   }
672: }

674: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
675: {
676:   PetscInt    ngp;
677:   PetscScalar gp_xi[GAUSS_POINTS][2];
678:   PetscScalar gp_weight[GAUSS_POINTS];
679:   PetscInt    p,i,j;
680:   PetscScalar Ni_p[NODES_PER_EL];
681:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
682:   PetscScalar J_p,fac,eta_avg;


685:   /* define quadrature rule */
686:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

688:   /* evaluate integral */
689:   for (p = 0; p < ngp; p++) {
690:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
691:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
692:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
693:     fac = gp_weight[p]*J_p;

695:     for (i = 0; i < NODES_PER_EL; i++) {
696:       for (j = 0; j < NODES_PER_EL; j++) {
697:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
698:       }
699:     }
700:   }

702:   /* scale */
703:   eta_avg = 0.0;
704:   for (p = 0; p < ngp; p++) {
705:     eta_avg += eta[p];
706:   }
707:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
708:   fac     = 1.0/eta_avg;
709:   for (i = 0; i < NODES_PER_EL; i++) {
710:     for (j = 0; j < NODES_PER_EL; j++) {
711:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
712:     }
713:   }
714: }

716: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
717: {
718:   PetscInt    ngp;
719:   PetscScalar gp_xi[GAUSS_POINTS][2];
720:   PetscScalar gp_weight[GAUSS_POINTS];
721:   PetscInt    p,i;
722:   PetscScalar Ni_p[NODES_PER_EL];
723:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
724:   PetscScalar J_p,fac;


727:   /* define quadrature rule */
728:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

730:   /* evaluate integral */
731:   for (p = 0; p < ngp; p++) {
732:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
733:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
734:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
735:     fac = gp_weight[p]*J_p;

737:     for (i = 0; i < NODES_PER_EL; i++) {
738:       Fe[NSD*i  ] += fac*Ni_p[i]*fx[p];
739:       Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
740:     }
741:   }
742: }

746: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
747: {
749:   /* get coords for the element */
750:   el_coords[NSD*0+0] = _coords[ej][ei].x;el_coords[NSD*0+1] = _coords[ej][ei].y;
751:   el_coords[NSD*1+0] = _coords[ej+1][ei].x;el_coords[NSD*1+1] = _coords[ej+1][ei].y;
752:   el_coords[NSD*2+0] = _coords[ej+1][ei+1].x;el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
753:   el_coords[NSD*3+0] = _coords[ej][ei+1].x;el_coords[NSD*3+1] = _coords[ej][ei+1].y;
754:   return(0);
755: }

759: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
760: {
761:   DM                     cda;
762:   Vec                    coords;
763:   DMDACoor2d               **_coords;
764:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
765:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
766:   PetscInt               sex,sey,mx,my;
767:   PetscInt               ei,ej;
768:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
769:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
770:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
771:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
772:   PetscScalar            el_coords[NODES_PER_EL*NSD];
773:   Vec                    local_properties;
774:   GaussPointCoefficients **props;
775:   PetscScalar            *prop_eta;
776:   PetscErrorCode         ierr;


780:   /* setup for coords */
781:   DMDAGetCoordinateDA(stokes_da,&cda);
782:   DMDAGetGhostedCoordinates(stokes_da,&coords);
783:   DMDAVecGetArray(cda,coords,&_coords);

785:   /* setup for coefficients */
786:   DMCreateLocalVector(properties_da,&local_properties);
787:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
788:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
789:   DMDAVecGetArray(properties_da,local_properties,&props);

791:   DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
792:   for (ej = sey; ej < sey+my; ej++) {
793:     for (ei = sex; ei < sex+mx; ei++) {
794:       /* get coords for the element */
795:       GetElementCoords(_coords,ei,ej,el_coords);

797:       /* get coefficients for the element */
798:       prop_eta = props[ej][ei].eta;

800:       /* initialise element stiffness matrix */
801:       PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
802:       PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
803:       PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
804:       PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

806:       /* form element stiffness matrix */
807:       FormStressOperatorQ1(Ae,el_coords,prop_eta);
808:       FormGradientOperatorQ1(Ge,el_coords);
809:       FormDivergenceOperatorQ1(De,el_coords);
810:       FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);

812:       /* insert element matrix into global matrix */
813:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
814:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
815:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
816:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
817:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
818:     }
819:   }
820:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
821:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

823:   DMDAVecRestoreArray(cda,coords,&_coords);

825:   DMDAVecRestoreArray(properties_da,local_properties,&props);
826:   VecDestroy(&local_properties);
827:   return(0);
828: }

832: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
833: {
834:   DM                     cda;
835:   Vec                    coords;
836:   DMDACoor2d               **_coords;
837:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
838:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
839:   PetscInt               sex,sey,mx,my;
840:   PetscInt               ei,ej;
841:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
842:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
843:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
844:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
845:   PetscScalar            el_coords[NODES_PER_EL*NSD];
846:   Vec                    local_properties;
847:   GaussPointCoefficients **props;
848:   PetscScalar            *prop_eta;
849:   PetscErrorCode         ierr;

852:   /* setup for coords */
853:   DMDAGetCoordinateDA(stokes_da,&cda);
854:   DMDAGetGhostedCoordinates(stokes_da,&coords);
855:   DMDAVecGetArray(cda,coords,&_coords);

857:   /* setup for coefficients */
858:   DMCreateLocalVector(properties_da,&local_properties);
859:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
860:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
861:   DMDAVecGetArray(properties_da,local_properties,&props);

863:   DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
864:   for (ej = sey; ej < sey+my; ej++) {
865:     for (ei = sex; ei < sex+mx; ei++) {
866:       /* get coords for the element */
867:       GetElementCoords(_coords,ei,ej,el_coords);

869:       /* get coefficients for the element */
870:       prop_eta = props[ej][ei].eta;

872:       /* initialise element stiffness matrix */
873:       PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
874:       PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
875:       PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
876:       PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);


879:       /* form element stiffness matrix */
880:       FormStressOperatorQ1(Ae,el_coords,prop_eta);
881:       FormGradientOperatorQ1(Ge,el_coords);
882:       /*               FormDivergenceOperatorQ1( De, el_coords ); */
883:       FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);

885:       /* insert element matrix into global matrix */
886:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
887:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
888:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
889:       /*     MatSetValuesStencil( A, NODES_PER_EL*P_DOFS,p_eqn, NODES_PER_EL*U_DOFS,u_eqn, De, ADD_VALUES ); */
890:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
891:     }
892:   }
893:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
894:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

896:   DMDAVecRestoreArray(cda,coords,&_coords);

898:   DMDAVecRestoreArray(properties_da,local_properties,&props);
899:   VecDestroy(&local_properties);
900:   return(0);
901: }

905: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
906: {
907:   PetscInt n;

910:   for (n = 0; n < 4; n++) {
911:     fields_F[ u_eqn[2*n  ].j ][ u_eqn[2*n  ].i ].u_dof = fields_F[ u_eqn[2*n  ].j ][ u_eqn[2*n  ].i ].u_dof+Fe_u[2*n  ];
912:     fields_F[ u_eqn[2*n+1].j ][ u_eqn[2*n+1].i ].v_dof = fields_F[ u_eqn[2*n+1].j ][ u_eqn[2*n+1].i ].v_dof+Fe_u[2*n+1];
913:     fields_F[ p_eqn[n].j     ][ p_eqn[n].i     ].p_dof = fields_F[ p_eqn[n].j     ][ p_eqn[n].i     ].p_dof+Fe_p[n    ];
914:   }
915:   return(0);
916: }

920: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
921: {
922:   DM                     cda;
923:   Vec                    coords;
924:   DMDACoor2d               **_coords;
925:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
926:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
927:   PetscInt               sex,sey,mx,my;
928:   PetscInt               ei,ej;
929:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
930:   PetscScalar            He[NODES_PER_EL*P_DOFS];
931:   PetscScalar            el_coords[NODES_PER_EL*NSD];
932:   Vec                    local_properties;
933:   GaussPointCoefficients **props;
934:   PetscScalar            *prop_fx,*prop_fy;
935:   Vec                    local_F;
936:   StokesDOF              **ff;
937:   PetscErrorCode         ierr;

940:   /* setup for coords */
941:   DMDAGetCoordinateDA(stokes_da,&cda);
942:   DMDAGetGhostedCoordinates(stokes_da,&coords);
943:   DMDAVecGetArray(cda,coords,&_coords);

945:   /* setup for coefficients */
946:   DMGetLocalVector(properties_da,&local_properties);
947:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
948:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
949:   DMDAVecGetArray(properties_da,local_properties,&props);

951:   /* get acces to the vector */
952:   DMGetLocalVector(stokes_da,&local_F);
953:   VecZeroEntries(local_F);
954:   DMDAVecGetArray(stokes_da,local_F,&ff);


957:   DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
958:   for (ej = sey; ej < sey+my; ej++) {
959:     for (ei = sex; ei < sex+mx; ei++) {
960:       /* get coords for the element */
961:       GetElementCoords(_coords,ei,ej,el_coords);

963:       /* get coefficients for the element */
964:       prop_fx = props[ej][ei].fx;
965:       prop_fy = props[ej][ei].fy;

967:       /* initialise element stiffness matrix */
968:       PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
969:       PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);


972:       /* form element stiffness matrix */
973:       FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);

975:       /* insert element matrix into global matrix */
976:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);

978:       DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
979:     }
980:   }

982:   DMDAVecRestoreArray(stokes_da,local_F,&ff);
983:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
984:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
985:   DMRestoreLocalVector(stokes_da,&local_F);


988:   DMDAVecRestoreArray(cda,coords,&_coords);

990:   DMDAVecRestoreArray(properties_da,local_properties,&props);
991:   DMRestoreLocalVector(properties_da,&local_properties);
992:   return(0);
993: }

997: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,
998:                                     PetscInt mx,PetscInt my,
999:                                     DM *_da,Vec *_X)
1000: {
1001:   DM             da,cda;
1002:   Vec            X,local_X;
1003:   StokesDOF      **_stokes;
1004:   Vec            coords;
1005:   DMDACoor2d       **_coords;
1006:   PetscInt       si,sj,ei,ej,i,j;

1010:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1011:                     mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,3,1,PETSC_NULL,PETSC_NULL,&da);
1012:   DMDASetFieldName(da,0,"anlytic_Vx");
1013:   DMDASetFieldName(da,1,"anlytic_Vy");
1014:   DMDASetFieldName(da,2,"analytic_P");


1017:   DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,PETSC_NULL,PETSC_NULL);


1020:   DMDAGetGhostedCoordinates(da,&coords);
1021:   DMDAGetCoordinateDA(da,&cda);
1022:   DMDAVecGetArray(cda,coords,&_coords);

1024:   DMCreateGlobalVector(da,&X);
1025:   DMCreateLocalVector(da,&local_X);
1026:   DMDAVecGetArray(da,local_X,&_stokes);

1028:   DMDAGetGhostCorners(da,&si,&sj,0,&ei,&ej,0);
1029:   for (j = sj; j < sj+ej; j++) {
1030:     for (i = si; i < si+ei; i++) {
1031:       double pos[2],pressure,vel[2],total_stress[3],strain_rate[3];

1033:       pos[0] = PetscRealPart(_coords[j][i].x);
1034:       pos[1] = PetscRealPart(_coords[j][i].y);

1036:       evaluate_solCx(pos,eta0,eta1,xc,nz,vel,&pressure,total_stress,strain_rate);

1038:       _stokes[j][i].u_dof = vel[0];
1039:       _stokes[j][i].v_dof = vel[1];
1040:       _stokes[j][i].p_dof = pressure;
1041:     }
1042:   }
1043:   DMDAVecRestoreArray(da,local_X,&_stokes);
1044:   DMDAVecRestoreArray(cda,coords,&_coords);

1046:   DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1047:   DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);

1049:   VecDestroy(&local_X);

1051:   *_da = da;
1052:   *_X  = X;

1054:   return(0);
1055: }

1059: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
1060: {
1062:   /* get the nodal fields */
1063:   nodal_fields[0].u_dof = fields[ej  ][ei  ].u_dof;nodal_fields[0].v_dof = fields[ej  ][ei  ].v_dof;nodal_fields[0].p_dof = fields[ej  ][ei  ].p_dof;
1064:   nodal_fields[1].u_dof = fields[ej+1][ei  ].u_dof;nodal_fields[1].v_dof = fields[ej+1][ei  ].v_dof;nodal_fields[1].p_dof = fields[ej+1][ei  ].p_dof;
1065:   nodal_fields[2].u_dof = fields[ej+1][ei+1].u_dof;nodal_fields[2].v_dof = fields[ej+1][ei+1].v_dof;nodal_fields[2].p_dof = fields[ej+1][ei+1].p_dof;
1066:   nodal_fields[3].u_dof = fields[ej  ][ei+1].u_dof;nodal_fields[3].v_dof = fields[ej  ][ei+1].v_dof;nodal_fields[3].p_dof = fields[ej  ][ei+1].p_dof;
1067:   return(0);
1068: }

1072: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
1073: {
1074:   DM          cda;
1075:   Vec         coords,X_analytic_local,X_local;
1076:   DMDACoor2d    **_coords;
1077:   PetscInt    sex,sey,mx,my;
1078:   PetscInt    ei,ej;
1079:   PetscScalar el_coords[NODES_PER_EL*NSD];
1080:   StokesDOF   **stokes_analytic,**stokes;
1081:   StokesDOF   stokes_analytic_e[4],stokes_e[4];

1083:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1084:   PetscScalar    Ni_p[NODES_PER_EL];
1085:   PetscInt       ngp;
1086:   PetscScalar    gp_xi[GAUSS_POINTS][2];
1087:   PetscScalar    gp_weight[GAUSS_POINTS];
1088:   PetscInt       p,i;
1089:   PetscScalar    J_p,fac;
1090:   PetscScalar    h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1091:   PetscInt       M;
1092:   PetscReal      xymin[2],xymax[2];

1096:   /* define quadrature rule */
1097:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

1099:   /* setup for coords */
1100:   DMDAGetCoordinateDA(stokes_da,&cda);
1101:   DMDAGetGhostedCoordinates(stokes_da,&coords);
1102:   DMDAVecGetArray(cda,coords,&_coords);

1104:   /* setup for analytic */
1105:   DMCreateLocalVector(stokes_da,&X_analytic_local);
1106:   DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1107:   DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1108:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1110:   /* setup for solution */
1111:   DMCreateLocalVector(stokes_da,&X_local);
1112:   DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1113:   DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1114:   DMDAVecGetArray(stokes_da,X_local,&stokes);

1116:   DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1117:   DMDAGetBoundingBox(stokes_da,xymin,xymax);

1119:   h = (xymax[0]-xymin[0])/((double)M);

1121:   tp_L2 = tu_L2 = tu_H1 = 0.0;

1123:   DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
1124:   for (ej = sey; ej < sey+my; ej++) {
1125:     for (ei = sex; ei < sex+mx; ei++) {
1126:       /* get coords for the element */
1127:       GetElementCoords(_coords,ei,ej,el_coords);
1128:       StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1129:       StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);

1131:       /* evaluate integral */
1132:       p_e_L2 = 0.0;
1133:       u_e_L2 = 0.0;
1134:       u_e_H1 = 0.0;
1135:       for (p = 0; p < ngp; p++) {
1136:         ConstructQ12D_Ni(gp_xi[p],Ni_p);
1137:         ConstructQ12D_GNi(gp_xi[p],GNi_p);
1138:         ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1139:         fac = gp_weight[p]*J_p;

1141:         for (i = 0; i < NODES_PER_EL; i++) {
1142:           PetscScalar u_error,v_error;

1144:           p_e_L2 = p_e_L2+fac*Ni_p[i]*(stokes_e[i].p_dof-stokes_analytic_e[i].p_dof)*(stokes_e[i].p_dof-stokes_analytic_e[i].p_dof);

1146:           u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1147:           v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1148:           u_e_L2  += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);

1150:           u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error              /* du/dx */
1151:                                +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error               /* du/dy */
1152:                                +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error               /* dv/dx */
1153:                                +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error);             /* dv/dy */
1154:         }
1155:       }

1157:       tp_L2 += p_e_L2;
1158:       tu_L2 += u_e_L2;
1159:       tu_H1 += u_e_H1;
1160:     }
1161:   }
1162:   MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1163:   MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1164:   MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1165:   p_L2 = PetscSqrtScalar(p_L2);
1166:   u_L2 = PetscSqrtScalar(u_L2);
1167:   u_H1 = PetscSqrtScalar(u_H1);

1169:   PetscPrintf(PETSC_COMM_WORLD,"%1.4e   %1.4e   %1.4e   %1.4e \n",PetscRealPart(h),PetscRealPart(p_L2),PetscRealPart(u_L2),PetscRealPart(u_H1));


1172:   DMDAVecRestoreArray(cda,coords,&_coords);

1174:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1175:   VecDestroy(&X_analytic_local);
1176:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1177:   VecDestroy(&X_local);
1178:   return(0);
1179: }

1183: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1184: {
1185:   DM                     da_Stokes,da_prop;
1186:   PetscInt               u_dof,p_dof,dof,stencil_width;
1187:   Mat                    A,B;
1188:   PetscInt               mxl,myl;
1189:   DM                     prop_cda,vel_cda;
1190:   Vec                    prop_coords,vel_coords;
1191:   PetscInt               si,sj,nx,ny,i,j,p;
1192:   Vec                    f,X;
1193:   PetscInt               prop_dof,prop_stencil_width;
1194:   Vec                    properties,l_properties;
1195:   PetscReal              dx,dy;
1196:   PetscInt               M,N;
1197:   DMDACoor2d               **_prop_coords,**_vel_coords;
1198:   GaussPointCoefficients **element_props;
1199:   PetscInt               its;
1200:   KSP                    ksp_S;
1201:   PetscInt               coefficient_structure = 0;
1202:   PetscInt               cpu_x,cpu_y,*lx = PETSC_NULL,*ly = PETSC_NULL;
1203:   PetscBool              use_gp_coords = PETSC_FALSE,set;
1204:   char                   filename[PETSC_MAX_PATH_LEN];
1205:   PetscErrorCode         ierr;

1208:   /* Generate the da for velocity and pressure */
1209:   /*
1210:   We use Q1 elements for the temperature.
1211:   FEM has a 9-point stencil (BOX) or connectivity pattern
1212:   Num nodes in each direction is mx+1, my+1
1213:   */
1214:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1215:   p_dof         = P_DOFS; /* p - pressure */
1216:   dof           = u_dof+p_dof;
1217:   stencil_width = 1;
1218:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1219:                              mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,PETSC_NULL,PETSC_NULL,&da_Stokes);
1220:   DMDASetFieldName(da_Stokes,0,"Vx");
1221:   DMDASetFieldName(da_Stokes,1,"Vy");
1222:   DMDASetFieldName(da_Stokes,2,"P");

1224:   /* unit box [0,1] x [0,1] */
1225:   DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,PETSC_NULL,PETSC_NULL);


1228:   /* Generate element properties, we will assume all material properties are constant over the element */
1229:   /* local number of elements */
1230:   DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,PETSC_NULL);

1232:   /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! // */
1233:   DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1234:   DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);

1236:   prop_dof           = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1237:   prop_stencil_width = 0;
1238:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1239:                                   mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
1240:   PetscFree(lx);
1241:   PetscFree(ly);

1243:   /* define centroid positions */
1244:   DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1245:   dx   = 1.0/((PetscReal)(M));
1246:   dy   = 1.0/((PetscReal)(N));

1248:   DMDASetUniformCoordinates(da_prop,0.0+0.5*dx,1.0-0.5*dx,0.0+0.5*dy,1.0-0.5*dy,PETSC_NULL,PETSC_NULL);

1250:   /* define coefficients */
1251:   PetscOptionsGetInt(PETSC_NULL,"-c_str",&coefficient_structure,PETSC_NULL);
1252:   /*     PetscPrintf( PETSC_COMM_WORLD, "Using coeficient structure %D \n", coefficient_structure ); */

1254:   DMCreateGlobalVector(da_prop,&properties);
1255:   DMCreateLocalVector(da_prop,&l_properties);
1256:   DMDAVecGetArray(da_prop,l_properties,&element_props);

1258:   DMDAGetCoordinateDA(da_prop,&prop_cda);
1259:   DMDAGetGhostedCoordinates(da_prop,&prop_coords);
1260:   DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);

1262:   DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);

1264:   DMDAGetCoordinateDA(da_Stokes,&vel_cda);
1265:   DMDAGetGhostedCoordinates(da_Stokes,&vel_coords);
1266:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);


1269:   /* interpolate the coordinates */
1270:   for (j = sj; j < sj+ny; j++) {
1271:     for (i = si; i < si+nx; i++) {
1272:       PetscInt    ngp;
1273:       PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1274:       PetscScalar el_coords[8];

1276:       GetElementCoords(_vel_coords,i,j,el_coords);
1277:       ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

1279:       for (p = 0; p < GAUSS_POINTS; p++) {
1280:         PetscScalar gp_x,gp_y;
1281:         PetscInt    n;
1282:         PetscScalar xi_p[2],Ni_p[4];

1284:         xi_p[0] = gp_xi[p][0];
1285:         xi_p[1] = gp_xi[p][1];
1286:         ConstructQ12D_Ni(xi_p,Ni_p);

1288:         gp_x = 0.0;
1289:         gp_y = 0.0;
1290:         for (n = 0; n < NODES_PER_EL; n++) {
1291:           gp_x = gp_x+Ni_p[n]*el_coords[2*n  ];
1292:           gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1293:         }
1294:         element_props[j][i].gp_coords[2*p  ] = gp_x;
1295:         element_props[j][i].gp_coords[2*p+1] = gp_y;
1296:       }
1297:     }
1298:   }

1300:   /* define the coefficients */
1301:   PetscOptionsGetBool(PETSC_NULL,"-use_gp_coords",&use_gp_coords,0);

1303:   for (j = sj; j < sj+ny; j++) {
1304:     for (i = si; i < si+nx; i++) {
1305:       PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1306:       PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1307:       PetscReal coord_x,coord_y;

1309:       if (coefficient_structure == 0) {
1310:         PetscReal opts_eta0,opts_eta1,opts_xc;
1311:         PetscInt  opts_nz;

1313:         opts_eta0 = 1.0;
1314:         opts_eta1 = 1.0;
1315:         opts_xc   = 0.5;
1316:         opts_nz   = 1;

1318:         PetscOptionsGetReal(PETSC_NULL,"-solcx_eta0",&opts_eta0,0);
1319:         PetscOptionsGetReal(PETSC_NULL,"-solcx_eta1",&opts_eta1,0);
1320:         PetscOptionsGetReal(PETSC_NULL,"-solcx_xc",&opts_xc,0);
1321:         PetscOptionsGetInt(PETSC_NULL,"-solcx_nz",&opts_nz,0);

1323:         for (p = 0; p < GAUSS_POINTS; p++) {
1324:           coord_x = centroid_x;
1325:           coord_y = centroid_y;
1326:           if (use_gp_coords == PETSC_TRUE) {
1327:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1328:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1329:           }


1332:           element_props[j][i].eta[p] = opts_eta0;
1333:           if (coord_x > opts_xc) {
1334:             element_props[j][i].eta[p] = opts_eta1;
1335:           }

1337:           element_props[j][i].fx[p] = 0.0;
1338:           element_props[j][i].fy[p] = sin((double)opts_nz*PETSC_PI*coord_y)*cos(1.0*PETSC_PI*coord_x);
1339:         }
1340:       } else if (coefficient_structure == 1) { /* square sinker */
1341:         PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;

1343:         opts_eta0 = 1.0;
1344:         opts_eta1 = 1.0;
1345:         opts_dx   = 0.50;
1346:         opts_dy   = 0.50;

1348:         PetscOptionsGetReal(PETSC_NULL,"-sinker_eta0",&opts_eta0,0);
1349:         PetscOptionsGetReal(PETSC_NULL,"-sinker_eta1",&opts_eta1,0);
1350:         PetscOptionsGetReal(PETSC_NULL,"-sinker_dx",&opts_dx,0);
1351:         PetscOptionsGetReal(PETSC_NULL,"-sinker_dy",&opts_dy,0);


1354:         for (p = 0; p < GAUSS_POINTS; p++) {
1355:           coord_x = centroid_x;
1356:           coord_y = centroid_y;
1357:           if (use_gp_coords == PETSC_TRUE) {
1358:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1359:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1360:           }

1362:           element_props[j][i].eta[p] = opts_eta0;
1363:           element_props[j][i].fx[p]  = 0.0;
1364:           element_props[j][i].fy[p]  = 0.0;

1366:           if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1367:             if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1368:               element_props[j][i].eta[p] =  opts_eta1;
1369:               element_props[j][i].fx[p]  =  0.0;
1370:               element_props[j][i].fy[p]  = -1.0;
1371:             }
1372:           }
1373:         }
1374:       } else if (coefficient_structure == 2) { /* circular sinker */
1375:         PetscReal opts_eta0,opts_eta1,opts_r,radius2;

1377:         opts_eta0 = 1.0;
1378:         opts_eta1 = 1.0;
1379:         opts_r    = 0.25;

1381:         PetscOptionsGetReal(PETSC_NULL,"-sinker_eta0",&opts_eta0,0);
1382:         PetscOptionsGetReal(PETSC_NULL,"-sinker_eta1",&opts_eta1,0);
1383:         PetscOptionsGetReal(PETSC_NULL,"-sinker_r",&opts_r,0);

1385:         for (p = 0; p < GAUSS_POINTS; p++) {
1386:           coord_x = centroid_x;
1387:           coord_y = centroid_y;
1388:           if (use_gp_coords == PETSC_TRUE) {
1389:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1390:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1391:           }

1393:           element_props[j][i].eta[p] = opts_eta0;
1394:           element_props[j][i].fx[p]  = 0.0;
1395:           element_props[j][i].fy[p]  = 0.0;

1397:           radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1398:           if (radius2 < opts_r*opts_r) {
1399:             element_props[j][i].eta[p] =  opts_eta1;
1400:             element_props[j][i].fx[p]  =  0.0;
1401:             element_props[j][i].fy[p]  = -1.0;
1402:           }
1403:         }
1404:       } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1405:         PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;

1407:         opts_eta0 = 1.0;
1408:         opts_eta1 = 1.0;
1409:         opts_r    = 0.25;
1410:         opts_c0x  = 0.35;       /* circle center */
1411:         opts_c0y  = 0.35;
1412:         opts_s0x  = 0.7;       /* square center */
1413:         opts_s0y  = 0.7;
1414:         opts_dx   = 0.25;
1415:         opts_dy   = 0.25;
1416:         opts_phi  = 25;

1418:         PetscOptionsGetReal(PETSC_NULL,"-sinker_eta0",&opts_eta0,0);
1419:         PetscOptionsGetReal(PETSC_NULL,"-sinker_eta1",&opts_eta1,0);
1420:         PetscOptionsGetReal(PETSC_NULL,"-sinker_r",&opts_r,0);
1421:         PetscOptionsGetReal(PETSC_NULL,"-sinker_c0x",&opts_c0x,0);
1422:         PetscOptionsGetReal(PETSC_NULL,"-sinker_c0y",&opts_c0y,0);
1423:         PetscOptionsGetReal(PETSC_NULL,"-sinker_s0x",&opts_s0x,0);
1424:         PetscOptionsGetReal(PETSC_NULL,"-sinker_s0y",&opts_s0y,0);
1425:         PetscOptionsGetReal(PETSC_NULL,"-sinker_dx",&opts_dx,0);
1426:         PetscOptionsGetReal(PETSC_NULL,"-sinker_dy",&opts_dy,0);
1427:         PetscOptionsGetReal(PETSC_NULL,"-sinker_phi",&opts_phi,0);
1428:         opts_phi *= PETSC_PI / 180;

1430:         for (p = 0; p < GAUSS_POINTS; p++) {
1431:           coord_x = centroid_x;
1432:           coord_y = centroid_y;
1433:           if (use_gp_coords == PETSC_TRUE) {
1434:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1435:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1436:           }

1438:           element_props[j][i].eta[p] = opts_eta0;
1439:           element_props[j][i].fx[p]  = 0.0;
1440:           element_props[j][i].fy[p]  = -0.2;

1442:           radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1443:           if (radius2 < opts_r*opts_r
1444:               || (   PetscAbs(+(coord_x - opts_s0x)*cos(opts_phi) + (coord_y - opts_s0y)*sin(opts_phi)) < opts_dx/2
1445:                   && PetscAbs(-(coord_x - opts_s0x)*sin(opts_phi) + (coord_y - opts_s0y)*cos(opts_phi)) < opts_dy/2)) {
1446:             element_props[j][i].eta[p] =  opts_eta1;
1447:             element_props[j][i].fx[p]  =  0.0;
1448:             element_props[j][i].fy[p]  = -1.0;
1449:           }
1450:         }
1451:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1452:     }
1453:   }
1454:   DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);

1456:   DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);

1458:   DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1459:   DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1460:   DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);


1463:   DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1464:   DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for Stokes eqn.","properties");


1467:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1468:   DMCreateMatrix(da_Stokes,MATAIJ,&A);
1469:   DMCreateMatrix(da_Stokes,MATAIJ,&B);
1470:   DMCreateGlobalVector(da_Stokes,&f);
1471:   DMCreateGlobalVector(da_Stokes,&X);

1473:   /* assemble A11 */
1474:   MatZeroEntries(A);
1475:   MatZeroEntries(B);
1476:   VecZeroEntries(f);

1478:   AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1479:   AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1480:   /* build force vector */
1481:   AssembleF_Stokes(f,da_Stokes,da_prop,properties);

1483:   DMDABCApplyFreeSlip(da_Stokes,A,f);
1484:   DMDABCApplyFreeSlip(da_Stokes,B,PETSC_NULL);

1486:   /* SOLVE */
1487:   KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1488:   KSPSetOptionsPrefix(ksp_S,"stokes_");
1489:   KSPSetOperators(ksp_S,A,B,SAME_NONZERO_PATTERN);
1490:   KSPSetDM(ksp_S,da_Stokes);
1491:   KSPSetDMActive(ksp_S,PETSC_FALSE);
1492:   KSPSetFromOptions(ksp_S);
1493:   {
1494:     PC pc;
1495:     const PetscInt ufields[] = {0,1},pfields[1] = {2};
1496:     KSPGetPC(ksp_S,&pc);
1497:     PCFieldSplitSetBlockSize(pc,3);
1498:     PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1499:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1500:   }

1502:   KSPSolve(ksp_S,f,X);

1504:   PetscOptionsGetString(PETSC_NULL,"-o",filename,sizeof filename,&set);
1505:   if (set) {
1506:     char        *ext;
1507:     PetscViewer viewer;
1508:     PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1509:     PetscStrrchr(filename,'.',&ext);
1510:     if (!strcmp("vts",ext)) {
1511:       PetscViewerSetType(viewer,PETSCVIEWERVTK);
1512:     } else {
1513:       PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1514:     }
1515:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1516:     PetscViewerFileSetName(viewer,filename);
1517:     VecView(X,viewer);
1518:     PetscViewerDestroy(&viewer);
1519:   }
1520:   DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");

1522:   KSPGetIterationNumber(ksp_S,&its);

1524:   /*
1525:   {
1526:   Vec x;
1527:   PetscInt L;
1528:   VecDuplicate(f,&x);
1529:   MatMult(A,X, x );
1530:   VecAXPY( x, -1.0, f );
1531:   VecNorm( x, NORM_2, &nrm );
1532:   PetscPrintf( PETSC_COMM_WORLD, "Its. %1.4d, norm |AX-f| = %1.5e \n", its, nrm );
1533:   VecDestroy(&x);

1535:   VecNorm( X, NORM_2, &nrm );
1536:   VecGetSize( X, &L );
1537:   PetscPrintf( PETSC_COMM_WORLD, "           norm |X|/sqrt{N} = %1.5e \n", nrm/sqrt( (PetscScalar)L ) );
1538:   }
1539:   */

1541:   if (coefficient_structure == 0) {
1542:     PetscReal   opts_eta0,opts_eta1,opts_xc;
1543:     PetscInt    opts_nz,N;
1544:     DM          da_Stokes_analytic;
1545:     Vec         X_analytic;
1546:     PetscReal   nrm1[3],nrm2[3],nrmI[3];

1548:     opts_eta0 = 1.0;
1549:     opts_eta1 = 1.0;
1550:     opts_xc   = 0.5;
1551:     opts_nz   = 1;

1553:     PetscOptionsGetReal(PETSC_NULL,"-solcx_eta0",&opts_eta0,0);
1554:     PetscOptionsGetReal(PETSC_NULL,"-solcx_eta1",&opts_eta1,0);
1555:     PetscOptionsGetReal(PETSC_NULL,"-solcx_xc",&opts_xc,0);
1556:     PetscOptionsGetInt(PETSC_NULL,"-solcx_nz",&opts_nz,0);


1559:     DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1560:     DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");

1562:     DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);


1565:     VecAXPY(X_analytic,-1.0,X);
1566:     VecGetSize(X_analytic,&N);
1567:     N    = N/3;

1569:     VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1570:     VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1571:     VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);

1573:     VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1574:     VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1575:     VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);

1577:     VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1578:     VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1579:     VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);

1581:     /*
1582:     PetscPrintf( PETSC_COMM_WORLD, "%1.4e   %1.4e %1.4e %1.4e    %1.4e %1.4e %1.4e    %1.4e %1.4e %1.4e\n",
1583:     1.0/mx,
1584:     nrm1[0]/(double)N, sqrt(nrm2[0]/(double)N), nrmI[0],
1585:     nrm1[1]/(double)N, sqrt(nrm2[1]/(double)N), nrmI[1],
1586:     nrm1[2]/(double)N, sqrt(nrm2[2]/(double)N), nrmI[2] );
1587:     */
1588:     DMDestroy(&da_Stokes_analytic);
1589:     VecDestroy(&X_analytic);
1590:   }


1593:   KSPDestroy(&ksp_S);
1594:   VecDestroy(&X);
1595:   VecDestroy(&f);
1596:   MatDestroy(&A);
1597:   MatDestroy(&B);

1599:   DMDestroy(&da_Stokes);
1600:   DMDestroy(&da_prop);

1602:   VecDestroy(&properties);
1603:   VecDestroy(&l_properties);

1605:   return(0);
1606: }

1610: int main(int argc,char **args)
1611: {
1613:   PetscInt       mx,my;

1615:   PetscInitialize(&argc,&args,(char *)0,help);

1617:   mx   = my = 10;
1618:   PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,PETSC_NULL);
1619:   PetscOptionsGetInt(PETSC_NULL,"-my",&my,PETSC_NULL);

1621:   solve_stokes_2d_coupled(mx,my);

1623:   PetscFinalize();
1624:   return 0;
1625: }

1627: /* -------------------------- helpers for boundary conditions -------------------------------- */

1631: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1632: {
1633:   DM             cda;
1634:   Vec            coords;
1635:   PetscInt       si,sj,nx,ny,i,j;
1636:   PetscInt       M,N;
1637:   DMDACoor2d       **_coords;
1638:   PetscInt       *g_idx;
1639:   PetscInt       *bc_global_ids;
1640:   PetscScalar    *bc_vals;
1641:   PetscInt       nbcs;
1642:   PetscInt       n_dofs;

1646:   /* enforce bc's */
1647:   DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);

1649:   DMDAGetCoordinateDA(da,&cda);
1650:   DMDAGetGhostedCoordinates(da,&coords);
1651:   DMDAVecGetArray(cda,coords,&_coords);
1652:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1653:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1655:   /* /// */

1657:   PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1658:   PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);

1660:   /* init the entries to -1 so VecSetValues will ignore them */
1661:   for (i = 0; i < ny*n_dofs; i++) {
1662:     bc_global_ids[i] = -1;
1663:   }

1665:   i = nx-1;
1666:   for (j = 0; j < ny; j++) {
1667:     PetscInt    local_id;

1669:     local_id = i+j*nx;

1671:     bc_global_ids[j] = g_idx[ n_dofs*local_id+d_idx ];

1673:     bc_vals[j] =  bc_val;
1674:   }
1675:   nbcs = 0;
1676:   if ((si+nx) == (M)) {
1677:     nbcs = ny;
1678:   }

1680:   if (b != PETSC_NULL) {
1681:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1682:     VecAssemblyBegin(b);
1683:     VecAssemblyEnd(b);
1684:   }
1685:   if (A != PETSC_NULL) {
1686:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1687:   }


1690:   PetscFree(bc_vals);
1691:   PetscFree(bc_global_ids);

1693:   DMDAVecRestoreArray(cda,coords,&_coords);
1694:   return(0);
1695: }

1699: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1700: {
1701:   DM             cda;
1702:   Vec            coords;
1703:   PetscInt       si,sj,nx,ny,i,j;
1704:   PetscInt       M,N;
1705:   DMDACoor2d       **_coords;
1706:   PetscInt       *g_idx;
1707:   PetscInt       *bc_global_ids;
1708:   PetscScalar    *bc_vals;
1709:   PetscInt       nbcs;
1710:   PetscInt       n_dofs;

1714:   /* enforce bc's */
1715:   DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);

1717:   DMDAGetCoordinateDA(da,&cda);
1718:   DMDAGetGhostedCoordinates(da,&coords);
1719:   DMDAVecGetArray(cda,coords,&_coords);
1720:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1721:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1723:   /* /// */

1725:   PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1726:   PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);

1728:   /* init the entries to -1 so VecSetValues will ignore them */
1729:   for (i = 0; i < ny*n_dofs; i++) {
1730:     bc_global_ids[i] = -1;
1731:   }

1733:   i = 0;
1734:   for (j = 0; j < ny; j++) {
1735:     PetscInt    local_id;

1737:     local_id = i+j*nx;

1739:     bc_global_ids[j] = g_idx[ n_dofs*local_id+d_idx ];

1741:     bc_vals[j] =  bc_val;
1742:   }
1743:   nbcs = 0;
1744:   if (si == 0) {
1745:     nbcs = ny;
1746:   }

1748:   if (b != PETSC_NULL) {
1749:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1750:     VecAssemblyBegin(b);
1751:     VecAssemblyEnd(b);
1752:   }
1753:   if (A != PETSC_NULL) {
1754:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1755:   }


1758:   PetscFree(bc_vals);
1759:   PetscFree(bc_global_ids);

1761:   DMDAVecRestoreArray(cda,coords,&_coords);
1762:   return(0);
1763: }

1767: static PetscErrorCode BCApply_NORTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1768: {
1769:   DM             cda;
1770:   Vec            coords;
1771:   PetscInt       si,sj,nx,ny,i,j;
1772:   PetscInt       M,N;
1773:   DMDACoor2d       **_coords;
1774:   PetscInt       *g_idx;
1775:   PetscInt       *bc_global_ids;
1776:   PetscScalar    *bc_vals;
1777:   PetscInt       nbcs;
1778:   PetscInt       n_dofs;

1782:   /* enforce bc's */
1783:   DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);

1785:   DMDAGetCoordinateDA(da,&cda);
1786:   DMDAGetGhostedCoordinates(da,&coords);
1787:   DMDAVecGetArray(cda,coords,&_coords);
1788:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1789:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1791:   /* /// */

1793:   PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1794:   PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);

1796:   /* init the entries to -1 so VecSetValues will ignore them */
1797:   for (i = 0; i < nx; i++) {
1798:     bc_global_ids[i] = -1;
1799:   }

1801:   j = ny-1;
1802:   for (i = 0; i < nx; i++) {
1803:     PetscInt    local_id;

1805:     local_id = i+j*nx;

1807:     bc_global_ids[i] = g_idx[ n_dofs*local_id+d_idx ];

1809:     bc_vals[i] =  bc_val;
1810:   }
1811:   nbcs = 0;
1812:   if ((sj+ny) == (N)) {
1813:     nbcs = nx;
1814:   }

1816:   if (b != PETSC_NULL) {
1817:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1818:     VecAssemblyBegin(b);
1819:     VecAssemblyEnd(b);
1820:   }
1821:   if (A != PETSC_NULL) {
1822:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1823:   }


1826:   PetscFree(bc_vals);
1827:   PetscFree(bc_global_ids);

1829:   DMDAVecRestoreArray(cda,coords,&_coords);
1830:   return(0);
1831: }

1835: static PetscErrorCode BCApply_SOUTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1836: {
1837:   DM             cda;
1838:   Vec            coords;
1839:   PetscInt       si,sj,nx,ny,i,j;
1840:   PetscInt       M,N;
1841:   DMDACoor2d       **_coords;
1842:   PetscInt       *g_idx;
1843:   PetscInt       *bc_global_ids;
1844:   PetscScalar    *bc_vals;
1845:   PetscInt       nbcs;
1846:   PetscInt       n_dofs;

1850:   /* enforce bc's */
1851:   DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);

1853:   DMDAGetCoordinateDA(da,&cda);
1854:   DMDAGetGhostedCoordinates(da,&coords);
1855:   DMDAVecGetArray(cda,coords,&_coords);
1856:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1857:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1859:   /* /// */

1861:   PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1862:   PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);

1864:   /* init the entries to -1 so VecSetValues will ignore them */
1865:   for (i = 0; i < nx; i++) {
1866:     bc_global_ids[i] = -1;
1867:   }

1869:   j = 0;
1870:   for (i = 0; i < nx; i++) {
1871:     PetscInt    local_id;

1873:     local_id = i+j*nx;

1875:     bc_global_ids[i] = g_idx[ n_dofs*local_id+d_idx ];

1877:     bc_vals[i] =  bc_val;
1878:   }
1879:   nbcs = 0;
1880:   if (sj == 0) {
1881:     nbcs = nx;
1882:   }


1885:   if (b != PETSC_NULL) {
1886:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1887:     VecAssemblyBegin(b);
1888:     VecAssemblyEnd(b);
1889:   }
1890:   if (A != PETSC_NULL) {
1891:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1892:   }


1895:   PetscFree(bc_vals);
1896:   PetscFree(bc_global_ids);

1898:   DMDAVecRestoreArray(cda,coords,&_coords);
1899:   return(0);
1900: }

1904: /*
1905: Free slip sides.
1906: */
1909: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1910: {

1914:   BCApply_NORTH(da_Stokes,1,0.0,A,f);
1915:   BCApply_EAST(da_Stokes,0,0.0,A,f);
1916:   BCApply_SOUTH(da_Stokes,1,0.0,A,f);
1917:   BCApply_WEST(da_Stokes,0,0.0,A,f);
1918:   return(0);
1919: }