Actual source code: ex43.c

petsc-master 2017-03-27
Report Typos and Errors
  1: static char help[] = "Solves the incompressible, variable viscosity Stokes equation in 2d on the unit domain \n\
  2: using Q1Q1 elements, stabilized with Bochev's polynomial projection method. \n\
  3: The models defined utilise free slip boundary conditions on all sides. \n\
  4: Options: \n"
  5: "\
  6:      -mx : Number of elements in the x-direction \n\
  7:      -my : Number of elements in the y-direction \n\
  8:      -o : Specify output filename for solution (will be petsc binary format or paraview format if the extension is .vts) \n\
  9:      -gnuplot : Output Gauss point coordinates, coefficients and u,p solution in gnuplot format \n\
 10:      -c_str : Indicates the structure of the coefficients to use \n"
 11: "\
 12:           -c_str 0 => Coefficient definition for an analytic solution with a vertical jump in viscosity at x = xc \n\
 13:                       This problem is driven by the forcing function f(x,y) = (0, sin(nz pi y)cos(pi x) \n\
 14:                          Parameters: \n\
 15:                               -solcx_eta0  : Viscosity to the left of the interface \n\
 16:                               -solcx_eta1  : Viscosity to the right of the interface \n\
 17:                               -solcx_xc    : Location of the interface \n\
 18:                               -solcx_nz    : Wavenumber in the y direction \n"
 19: "\
 20:           -c_str 1 => Coefficient definition for a dense rectangular blob located at the center of the domain \n\
 21:                          Parameters: \n\
 22:                               -sinker_eta0 : Viscosity of the background fluid \n\
 23:                               -sinker_eta1 : Viscosity of the blob \n\
 24:                               -sinker_dx   : Width of the blob \n\
 25:                               -sinker_dy   : Height of the blob \n"
 26: "\
 27:           -c_str 2 => Coefficient definition for a dense circular blob located at the center of the domain \n\
 28:                          Parameters: \n\
 29:                               -sinker_eta0 : Viscosity of the background fluid \n\
 30:                               -sinker_eta1 : Viscosity of the blob \n\
 31:                               -sinker_r    : Radius of the blob \n"
 32: "\
 33:           -c_str 3 => Coefficient definition for a dense circular and rectangular inclusion (located at the center of the domain) \n\
 34:                               -sinker_eta0 : Viscosity of the background fluid \n\
 35:                               -sinker_eta1 : Viscosity of the two inclusions \n\
 36:                               -sinker_r    : Radius of the circular inclusion \n\
 37:                               -sinker_c0x  : Origin (x-coord) of the circular inclusion \n\
 38:                               -sinker_c0y  : Origin (y-coord) of the circular inclusion \n\
 39:                               -sinker_dx   : Width of the rectangular inclusion \n\
 40:                               -sinker_dy   : Height of the rectangular inclusion \n\
 41:                               -sinker_phi  : Rotation angle of the rectangular inclusion \n\
 42:      -use_gp_coords : Evaluate the viscosity and force term at the global coordinates of each quadrature point \n\
 43:                       By default, the viscosity and force term are evaulated at the element center and applied as a constant over the entire element \n";

 45: /* Contributed by Dave May */

 47:  #include <petscksp.h>
 48:  #include <petscdm.h>
 49:  #include <petscdmda.h>

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

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


 57: #define NSD            2 /* number of spatial dimensions */
 58: #define NODES_PER_EL   4 /* nodes per element */
 59: #define U_DOFS         2 /* degrees of freedom per velocity node */
 60: #define P_DOFS         1 /* degrees of freedom per pressure node */
 61: #define GAUSS_POINTS   4

 63: /* Gauss point based evaluation 8+4+4+4 = 20 */
 64: typedef struct {
 65:   PetscScalar gp_coords[2*GAUSS_POINTS];
 66:   PetscScalar eta[GAUSS_POINTS];
 67:   PetscScalar fx[GAUSS_POINTS];
 68:   PetscScalar fy[GAUSS_POINTS];
 69: } GaussPointCoefficients;

 71: typedef struct {
 72:   PetscScalar u_dof;
 73:   PetscScalar v_dof;
 74:   PetscScalar p_dof;
 75: } StokesDOF;


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

 90:   Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
 91:   Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
 92:   Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
 93:   Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
 94: }

 96: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
 97: {
 98:   PetscScalar xi  = _xi[0];
 99:   PetscScalar eta = _xi[1];

101:   GNi[0][0] = -0.25*(1.0-eta);
102:   GNi[0][1] = -0.25*(1.0+eta);
103:   GNi[0][2] =   0.25*(1.0+eta);
104:   GNi[0][3] =   0.25*(1.0-eta);

106:   GNi[1][0] = -0.25*(1.0-xi);
107:   GNi[1][1] =   0.25*(1.0-xi);
108:   GNi[1][2] =   0.25*(1.0+xi);
109:   GNi[1][3] = -0.25*(1.0+xi);
110: }

112: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
113: {
114:   PetscScalar J00,J01,J10,J11,J;
115:   PetscScalar iJ00,iJ01,iJ10,iJ11;
116:   PetscInt    i;

118:   J00 = J01 = J10 = J11 = 0.0;
119:   for (i = 0; i < NODES_PER_EL; i++) {
120:     PetscScalar cx = coords[2*i];
121:     PetscScalar cy = coords[2*i+1];

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

130:   iJ00 =  J11/J;
131:   iJ01 = -J01/J;
132:   iJ10 = -J10/J;
133:   iJ11 =  J00/J;

135:   for (i = 0; i < NODES_PER_EL; i++) {
136:     GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
137:     GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
138:   }

140:   if (det_J) *det_J = J;
141: }

143: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
144: {
145:   *ngp         = 4;
146:   gp_xi[0][0]  = -0.57735026919; gp_xi[0][1] = -0.57735026919;
147:   gp_xi[1][0]  = -0.57735026919; gp_xi[1][1] =  0.57735026919;
148:   gp_xi[2][0]  =  0.57735026919; gp_xi[2][1] =  0.57735026919;
149:   gp_xi[3][0]  =  0.57735026919; gp_xi[3][1] = -0.57735026919;
150:   gp_weight[0] = 1.0;
151:   gp_weight[1] = 1.0;
152:   gp_weight[2] = 1.0;
153:   gp_weight[3] = 1.0;
154: }

156: /* procs to the left claim the ghost node as their element */
157: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
158: {
159:   PetscInt m,n,p,M,N,P;
160:   PetscInt sx,sy,sz;

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

166:   if (mxl) {
167:     *mxl = m;
168:     if ((sx+m) == M) *mxl = m-1;  /* last proc */
169:   }
170:   if (myl) {
171:     *myl = n;
172:     if ((sy+n) == N) *myl = n-1;  /* last proc */
173:   }
174:   if (mzl) {
175:     *mzl = p;
176:     if ((sz+p) == P) *mzl = p-1;  /* last proc */
177:   }
178:   return(0);
179: }

181: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
182: {
183:   PetscInt       si,sj,sk;

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

189:   *sx = si;
190:   if (si) *sx = si+1;

192:   *sy = sj;
193:   if (sj) *sy = sj+1;

195:   if (sz) {
196:     *sz = sk;
197:     if (sk) *sz = sk+1;
198:   }

200:   DMDAGetLocalElementSize(da,mx,my,mz);
201:   return(0);
202: }

204: /*
205: i,j are the element indices
206: The unknown is a vector quantity.
207: The s[].c is used to indicate the degree of freedom.
208: */
209: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
210: {
212:   /* velocity */
213:   /* node 0 */
214:   s_u[0].i = i; s_u[0].j = j; s_u[0].c = 0;       /* Vx0 */
215:   s_u[1].i = i; s_u[1].j = j; s_u[1].c = 1;       /* Vy0 */

217:   /* node 1 */
218:   s_u[2].i = i; s_u[2].j = j+1; s_u[2].c = 0;     /* Vx1 */
219:   s_u[3].i = i; s_u[3].j = j+1; s_u[3].c = 1;     /* Vy1 */

221:   /* node 2 */
222:   s_u[4].i = i+1; s_u[4].j = j+1; s_u[4].c = 0;   /* Vx2 */
223:   s_u[5].i = i+1; s_u[5].j = j+1; s_u[5].c = 1;   /* Vy2 */

225:   /* node 3 */
226:   s_u[6].i = i+1; s_u[6].j = j; s_u[6].c = 0;     /* Vx3 */
227:   s_u[7].i = i+1; s_u[7].j = j; s_u[7].c = 1;     /* Vy3 */

229:   /* pressure */
230:   s_p[0].i = i;   s_p[0].j = j;   s_p[0].c = 2; /* P0 */
231:   s_p[1].i = i;   s_p[1].j = j+1; s_p[1].c = 2; /* P1 */
232:   s_p[2].i = i+1; s_p[2].j = j+1; s_p[2].c = 2; /* P2 */
233:   s_p[3].i = i+1; s_p[3].j = j;   s_p[3].c = 2; /* P3 */
234:   return(0);
235: }

237: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
238: {
240:   PetscMPIInt    rank;
241:   PetscInt       proc_I,proc_J;
242:   PetscInt       cpu_x,cpu_y;
243:   PetscInt       local_mx,local_my;
244:   Vec            vlx,vly;
245:   PetscInt       *LX,*LY,i;
246:   PetscScalar    *_a;
247:   Vec            V_SEQ;
248:   VecScatter     ctx;

251:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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

255:   proc_J = rank/cpu_x;
256:   proc_I = rank-cpu_x*proc_J;

258:   PetscMalloc1(cpu_x,&LX);
259:   PetscMalloc1(cpu_y,&LY);

261:   DMDAGetLocalElementSize(da,&local_mx,&local_my,NULL);
262:   VecCreate(PETSC_COMM_WORLD,&vlx);
263:   VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
264:   VecSetFromOptions(vlx);

266:   VecCreate(PETSC_COMM_WORLD,&vly);
267:   VecSetSizes(vly,PETSC_DECIDE,cpu_y);
268:   VecSetFromOptions(vly);

270:   VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
271:   VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
272:   VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
273:   VecAssemblyBegin(vly);VecAssemblyEnd(vly);

275:   VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
276:   VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
277:   VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
278:   VecGetArray(V_SEQ,&_a);
279:   for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
280:   VecRestoreArray(V_SEQ,&_a);
281:   VecScatterDestroy(&ctx);
282:   VecDestroy(&V_SEQ);

284:   VecScatterCreateToAll(vly,&ctx,&V_SEQ);
285:   VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
286:   VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
287:   VecGetArray(V_SEQ,&_a);
288:   for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
289:   VecRestoreArray(V_SEQ,&_a);
290:   VecScatterDestroy(&ctx);
291:   VecDestroy(&V_SEQ);

293:   *_lx = LX;
294:   *_ly = LY;

296:   VecDestroy(&vlx);
297:   VecDestroy(&vly);
298:   return(0);
299: }

301: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
302: {
303:   DM             cda;
304:   Vec            coords;
305:   DMDACoor2d     **_coords;
306:   PetscInt       si,sj,nx,ny,i,j;
307:   FILE           *fp;
308:   char           fname[PETSC_MAX_PATH_LEN];
309:   PetscMPIInt    rank;

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

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

320:   DMGetCoordinateDM(da,&cda);
321:   DMGetCoordinatesLocal(da,&coords);
322:   DMDAVecGetArrayRead(cda,coords,&_coords);
323:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
324:   for (j = sj; j < sj+ny-1; j++) {
325:     for (i = si; i < si+nx-1; i++) {
326:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
327:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i].x),(double)PetscRealPart(_coords[j+1][i].y));
328:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i+1].x),(double)PetscRealPart(_coords[j+1][i+1].y));
329:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i+1].x),(double)PetscRealPart(_coords[j][i+1].y));
330:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
331:     }
332:   }
333:   DMDAVecRestoreArrayRead(cda,coords,&_coords);

335:   PetscFClose(PETSC_COMM_SELF,fp);
336:   return(0);
337: }

339: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
340: {
341:   DM             cda;
342:   Vec            coords,local_fields;
343:   DMDACoor2d     **_coords;
344:   FILE           *fp;
345:   char           fname[PETSC_MAX_PATH_LEN];
346:   PetscMPIInt    rank;
347:   PetscInt       si,sj,nx,ny,i,j;
348:   PetscInt       n_dofs,d;
349:   PetscScalar    *_fields;

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

358:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
359:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
360:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
361:   for (d = 0; d < n_dofs; d++) {
362:     const char *field_name;
363:     DMDAGetFieldName(da,d,&field_name);
364:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
365:   }
366:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");

368:   DMGetCoordinateDM(da,&cda);
369:   DMGetCoordinatesLocal(da,&coords);
370:   DMDAVecGetArray(cda,coords,&_coords);
371:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

373:   DMCreateLocalVector(da,&local_fields);
374:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
375:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
376:   VecGetArray(local_fields,&_fields);

378:   for (j = sj; j < sj+ny; j++) {
379:     for (i = si; i < si+nx; i++) {
380:       PetscScalar coord_x,coord_y;
381:       PetscScalar field_d;

383:       coord_x = _coords[j][i].x;
384:       coord_y = _coords[j][i].y;

386:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
387:       for (d = 0; d < n_dofs; d++) {
388:         field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
389:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
390:       }
391:       PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
392:     }
393:   }
394:   VecRestoreArray(local_fields,&_fields);
395:   VecDestroy(&local_fields);

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

399:   PetscFClose(PETSC_COMM_SELF,fp);
400:   return(0);
401: }

403: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
404: {
405:   DM                     cda;
406:   Vec                    local_fields;
407:   FILE                   *fp;
408:   char                   fname[PETSC_MAX_PATH_LEN];
409:   PetscMPIInt            rank;
410:   PetscInt               si,sj,nx,ny,i,j,p;
411:   PetscInt               n_dofs,d;
412:   GaussPointCoefficients **_coefficients;
413:   PetscErrorCode         ierr;

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

421:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
422:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
423:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
424:   for (d = 0; d < n_dofs; d++) {
425:     const char *field_name;
426:     DMDAGetFieldName(da,d,&field_name);
427:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
428:   }
429:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");

431:   DMGetCoordinateDM(da,&cda);
432:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

434:   DMCreateLocalVector(da,&local_fields);
435:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
436:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
437:   DMDAVecGetArray(da,local_fields,&_coefficients);

439:   for (j = sj; j < sj+ny; j++) {
440:     for (i = si; i < si+nx; i++) {
441:       PetscScalar coord_x,coord_y;

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

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

449:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e",(double)PetscRealPart(_coefficients[j][i].eta[p]),(double)PetscRealPart(_coefficients[j][i].fx[p]),(double)PetscRealPart(_coefficients[j][i].fy[p]));
450:         PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
451:       }
452:     }
453:   }
454:   DMDAVecRestoreArray(da,local_fields,&_coefficients);
455:   VecDestroy(&local_fields);

457:   PetscFClose(PETSC_COMM_SELF,fp);
458:   return(0);
459: }

461: static PetscInt ASS_MAP_wIwDI_uJuDJ(PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
462: {
463:   PetscInt ij;
464:   PetscInt r,c,nc;

466:   nc = u_NPE*u_dof;
467:   r = w_dof*wi+wd;
468:   c = u_dof*ui+ud;
469:   ij = r*nc+c;
470:   return ij;
471: }

473: /*
474:  D = [ 2.eta   0   0   ]
475:      [   0   2.eta 0   ]
476:      [   0     0   eta ]
477:  
478:  B = [ d_dx   0   ]
479:      [  0    d_dy ]
480:      [ d_dy  d_dx ]
481: */
482: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
483: {
484:   PetscInt       ngp;
485:   PetscScalar    gp_xi[GAUSS_POINTS][2];
486:   PetscScalar    gp_weight[GAUSS_POINTS];
487:   PetscInt       p,i,j,k;
488:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
489:   PetscScalar    J_p,tildeD[3];
490:   PetscScalar    B[3][U_DOFS*NODES_PER_EL];

492:   /* define quadrature rule */
493:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

495:   /* evaluate integral */
496:   for (p = 0; p < ngp; p++) {
497:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
498:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);

500:     for (i = 0; i < NODES_PER_EL; i++) {
501:       PetscScalar d_dx_i = GNx_p[0][i];
502:       PetscScalar d_dy_i = GNx_p[1][i];

504:       B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
505:       B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
506:       B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
507:     }

509:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
510:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
511:     tildeD[2] =       gp_weight[p]*J_p*eta[p];

513:     /* form Bt tildeD B */
514:     /*
515:     Ke_ij = Bt_ik . D_kl . B_lj
516:           = B_ki . D_kl . B_lj
517:           = B_ki . D_kk . B_kj
518:     */
519:     for (i = 0; i < 8; i++) {
520:       for (j = 0; j < 8; j++) {
521:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
522:           Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
523:         }
524:       }
525:     }
526:   }
527: }

529: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
530: {
531:   PetscInt    ngp;
532:   PetscScalar gp_xi[GAUSS_POINTS][2];
533:   PetscScalar gp_weight[GAUSS_POINTS];
534:   PetscInt    p,i,j,di;
535:   PetscScalar Ni_p[NODES_PER_EL];
536:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
537:   PetscScalar J_p,fac;

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

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

549:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
550:       for (di = 0; di < NSD; di++) { /* u dofs */
551:         for (j = 0; j < 4; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
552:           PetscInt IJ;
553:           IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);

555:           Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
556:         }
557:       }
558:     }
559:   }
560: }

562: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
563: {
564:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
565:   PetscInt    i,j;
566:   PetscInt    nr_g,nc_g;

568:   PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
569:   FormGradientOperatorQ1(Ge,coords);

571:   nr_g = U_DOFS*NODES_PER_EL;
572:   nc_g = P_DOFS*NODES_PER_EL;

574:   for (i = 0; i < nr_g; i++) {
575:     for (j = 0; j < nc_g; j++) {
576:       De[nr_g*j+i] = Ge[nc_g*i+j];
577:     }
578:   }
579: }

581: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
582: {
583:   PetscInt    ngp;
584:   PetscScalar gp_xi[GAUSS_POINTS][2];
585:   PetscScalar gp_weight[GAUSS_POINTS];
586:   PetscInt    p,i,j;
587:   PetscScalar Ni_p[NODES_PER_EL];
588:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
589:   PetscScalar J_p,fac,eta_avg;

591:   /* define quadrature rule */
592:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

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

601:     for (i = 0; i < NODES_PER_EL; i++) {
602:       for (j = 0; j < NODES_PER_EL; j++) {
603:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
604:       }
605:     }
606:   }

608:   /* scale */
609:   eta_avg = 0.0;
610:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
611:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
612:   fac     = 1.0/eta_avg;
613:   for (i = 0; i < NODES_PER_EL; i++) {
614:     for (j = 0; j < NODES_PER_EL; j++) {
615:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
616:     }
617:   }
618: }

620: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
621: {
622:   PetscInt    ngp;
623:   PetscScalar gp_xi[GAUSS_POINTS][2];
624:   PetscScalar gp_weight[GAUSS_POINTS];
625:   PetscInt    p,i,j;
626:   PetscScalar Ni_p[NODES_PER_EL];
627:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
628:   PetscScalar J_p,fac,eta_avg;

630:   /* define quadrature rule */
631:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

633:   /* evaluate integral */
634:   for (p = 0; p < ngp; p++) {
635:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
636:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
637:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
638:     fac = gp_weight[p]*J_p;

640:     for (i = 0; i < NODES_PER_EL; i++) {
641:       for (j = 0; j < NODES_PER_EL; j++) {
642:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
643:       }
644:     }
645:   }

647:   /* scale */
648:   eta_avg = 0.0;
649:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
650:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
651:   fac     = 1.0/eta_avg;
652:   for (i = 0; i < NODES_PER_EL; i++) {
653:     for (j = 0; j < NODES_PER_EL; j++) {
654:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
655:     }
656:   }
657: }

659: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
660: {
661:   PetscInt    ngp;
662:   PetscScalar gp_xi[GAUSS_POINTS][2];
663:   PetscScalar gp_weight[GAUSS_POINTS];
664:   PetscInt    p,i;
665:   PetscScalar Ni_p[NODES_PER_EL];
666:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
667:   PetscScalar J_p,fac;

669:   /* define quadrature rule */
670:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

672:   /* evaluate integral */
673:   for (p = 0; p < ngp; p++) {
674:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
675:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
676:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
677:     fac = gp_weight[p]*J_p;

679:     for (i = 0; i < NODES_PER_EL; i++) {
680:       Fe[NSD*i]   += fac*Ni_p[i]*fx[p];
681:       Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
682:     }
683:   }
684: }

686: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
687: {
689:   /* get coords for the element */
690:   el_coords[NSD*0] = _coords[ej][ei].x;     el_coords[NSD*0+1] = _coords[ej][ei].y;
691:   el_coords[NSD*1] = _coords[ej+1][ei].x;   el_coords[NSD*1+1] = _coords[ej+1][ei].y;
692:   el_coords[NSD*2] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
693:   el_coords[NSD*3] = _coords[ej][ei+1].x;   el_coords[NSD*3+1] = _coords[ej][ei+1].y;
694:   return(0);
695: }

697: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
698: {
699:   DM                     cda;
700:   Vec                    coords;
701:   DMDACoor2d             **_coords;
702:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
703:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
704:   PetscInt               sex,sey,mx,my;
705:   PetscInt               ei,ej;
706:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
707:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
708:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
709:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
710:   PetscScalar            el_coords[NODES_PER_EL*NSD];
711:   Vec                    local_properties;
712:   GaussPointCoefficients **props;
713:   PetscScalar            *prop_eta;
714:   PetscErrorCode         ierr;

717:   /* setup for coords */
718:   DMGetCoordinateDM(stokes_da,&cda);
719:   DMGetCoordinatesLocal(stokes_da,&coords);
720:   DMDAVecGetArray(cda,coords,&_coords);

722:   /* setup for coefficients */
723:   DMCreateLocalVector(properties_da,&local_properties);
724:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
725:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
726:   DMDAVecGetArray(properties_da,local_properties,&props);

728:   DMDAGetElementCorners(stokes_da,&sex,&sey,NULL,&mx,&my,NULL);
729:   for (ej = sey; ej < sey+my; ej++) {
730:     for (ei = sex; ei < sex+mx; ei++) {
731:       /* get coords for the element */
732:       GetElementCoords(_coords,ei,ej,el_coords);

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

737:       /* initialise element stiffness matrix */
738:       PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
739:       PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
740:       PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
741:       PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

743:       /* form element stiffness matrix */
744:       FormStressOperatorQ1(Ae,el_coords,prop_eta);
745:       FormGradientOperatorQ1(Ge,el_coords);
746:       FormDivergenceOperatorQ1(De,el_coords);
747:       FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);

749:       /* insert element matrix into global matrix */
750:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
751:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
752:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
753:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
754:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
755:     }
756:   }
757:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
758:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

762:   DMDAVecRestoreArray(properties_da,local_properties,&props);
763:   VecDestroy(&local_properties);
764:   return(0);
765: }

767: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
768: {
769:   DM                     cda;
770:   Vec                    coords;
771:   DMDACoor2d             **_coords;
772:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
773:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
774:   PetscInt               sex,sey,mx,my;
775:   PetscInt               ei,ej;
776:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
777:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
778:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
779:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
780:   PetscScalar            el_coords[NODES_PER_EL*NSD];
781:   Vec                    local_properties;
782:   GaussPointCoefficients **props;
783:   PetscScalar            *prop_eta;
784:   PetscErrorCode         ierr;

787:   /* setup for coords */
788:   DMGetCoordinateDM(stokes_da,&cda);
789:   DMGetCoordinatesLocal(stokes_da,&coords);
790:   DMDAVecGetArray(cda,coords,&_coords);

792:   /* setup for coefficients */
793:   DMCreateLocalVector(properties_da,&local_properties);
794:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
795:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
796:   DMDAVecGetArray(properties_da,local_properties,&props);

798:   DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
799:   for (ej = sey; ej < sey+my; ej++) {
800:     for (ei = sex; ei < sex+mx; ei++) {
801:       /* get coords for the element */
802:       GetElementCoords(_coords,ei,ej,el_coords);

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

807:       /* initialise element stiffness matrix */
808:       PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
809:       PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
810:       PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
811:       PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);

813:       /* form element stiffness matrix */
814:       FormStressOperatorQ1(Ae,el_coords,prop_eta);
815:       FormGradientOperatorQ1(Ge,el_coords);
816:       FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);

818:       /* insert element matrix into global matrix */
819:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
820:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
821:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
822:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
823:     }
824:   }
825:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
826:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

830:   DMDAVecRestoreArray(properties_da,local_properties,&props);
831:   VecDestroy(&local_properties);
832:   return(0);
833: }

835: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
836: {
837:   PetscInt n;

840:   for (n = 0; n < 4; n++) {
841:     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];
842:     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];
843:     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];
844:   }
845:   return(0);
846: }

848: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
849: {
850:   DM                     cda;
851:   Vec                    coords;
852:   DMDACoor2d             **_coords;
853:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
854:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
855:   PetscInt               sex,sey,mx,my;
856:   PetscInt               ei,ej;
857:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
858:   PetscScalar            He[NODES_PER_EL*P_DOFS];
859:   PetscScalar            el_coords[NODES_PER_EL*NSD];
860:   Vec                    local_properties;
861:   GaussPointCoefficients **props;
862:   PetscScalar            *prop_fx,*prop_fy;
863:   Vec                    local_F;
864:   StokesDOF              **ff;
865:   PetscErrorCode         ierr;

868:   /* setup for coords */
869:   DMGetCoordinateDM(stokes_da,&cda);
870:   DMGetCoordinatesLocal(stokes_da,&coords);
871:   DMDAVecGetArray(cda,coords,&_coords);

873:   /* setup for coefficients */
874:   DMGetLocalVector(properties_da,&local_properties);
875:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
876:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
877:   DMDAVecGetArray(properties_da,local_properties,&props);

879:   /* get acces to the vector */
880:   DMGetLocalVector(stokes_da,&local_F);
881:   VecZeroEntries(local_F);
882:   DMDAVecGetArray(stokes_da,local_F,&ff);

884:   DMDAGetElementCorners(stokes_da,&sex,&sey,NULL,&mx,&my,NULL);
885:   for (ej = sey; ej < sey+my; ej++) {
886:     for (ei = sex; ei < sex+mx; ei++) {
887:       /* get coords for the element */
888:       GetElementCoords(_coords,ei,ej,el_coords);

890:       /* get coefficients for the element */
891:       prop_fx = props[ej][ei].fx;
892:       prop_fy = props[ej][ei].fy;

894:       /* initialise element stiffness matrix */
895:       PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
896:       PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);

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

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

904:       DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
905:     }
906:   }

908:   DMDAVecRestoreArray(stokes_da,local_F,&ff);
909:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
910:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
911:   DMRestoreLocalVector(stokes_da,&local_F);

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

915:   DMDAVecRestoreArray(properties_da,local_properties,&props);
916:   DMRestoreLocalVector(properties_da,&local_properties);
917:   return(0);
918: }

920: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM *_da,Vec *_X)
921: {
922:   DM             da,cda;
923:   Vec            X;
924:   StokesDOF      **_stokes;
925:   Vec            coords;
926:   DMDACoor2d     **_coords;
927:   PetscInt       si,sj,ei,ej,i,j;

931:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,&da);
932:   DMSetFromOptions(da);
933:   DMSetUp(da);
934:   DMDASetFieldName(da,0,"anlytic_Vx");
935:   DMDASetFieldName(da,1,"anlytic_Vy");
936:   DMDASetFieldName(da,2,"analytic_P");

938:   DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.,0.);

940:   DMGetCoordinatesLocal(da,&coords);
941:   DMGetCoordinateDM(da,&cda);
942:   DMDAVecGetArray(cda,coords,&_coords);

944:   DMCreateGlobalVector(da,&X);
945:   DMDAVecGetArray(da,X,&_stokes);

947:   DMDAGetCorners(da,&si,&sj,0,&ei,&ej,0);
948:   for (j = sj; j < sj+ej; j++) {
949:     for (i = si; i < si+ei; i++) {
950:       PetscReal pos[2],pressure,vel[2],total_stress[3],strain_rate[3];

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

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

957:       _stokes[j][i].u_dof = vel[0];
958:       _stokes[j][i].v_dof = vel[1];
959:       _stokes[j][i].p_dof = pressure;
960:     }
961:   }
962:   DMDAVecRestoreArray(da,X,&_stokes);
963:   DMDAVecRestoreArray(cda,coords,&_coords);

965:   *_da = da;
966:   *_X  = X;
967:   return(0);
968: }

970: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
971: {
973:   /* get the nodal fields */
974:   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;
975:   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;
976:   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;
977:   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;
978:   return(0);
979: }

981: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
982: {
983:   DM             cda;
984:   Vec            coords,X_analytic_local,X_local;
985:   DMDACoor2d     **_coords;
986:   PetscInt       sex,sey,mx,my;
987:   PetscInt       ei,ej;
988:   PetscScalar    el_coords[NODES_PER_EL*NSD];
989:   StokesDOF      **stokes_analytic,**stokes;
990:   StokesDOF      stokes_analytic_e[4],stokes_e[4];

992:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
993:   PetscScalar    Ni_p[NODES_PER_EL];
994:   PetscInt       ngp;
995:   PetscScalar    gp_xi[GAUSS_POINTS][2];
996:   PetscScalar    gp_weight[GAUSS_POINTS];
997:   PetscInt       p,i;
998:   PetscScalar    J_p,fac;
999:   PetscScalar    h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1000:   PetscInt       M;
1001:   PetscReal      xymin[2],xymax[2];

1005:   /* define quadrature rule */
1006:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

1008:   /* setup for coords */
1009:   DMGetCoordinateDM(stokes_da,&cda);
1010:   DMGetCoordinatesLocal(stokes_da,&coords);
1011:   DMDAVecGetArray(cda,coords,&_coords);

1013:   /* setup for analytic */
1014:   DMCreateLocalVector(stokes_da,&X_analytic_local);
1015:   DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1016:   DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1017:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1019:   /* setup for solution */
1020:   DMCreateLocalVector(stokes_da,&X_local);
1021:   DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1022:   DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1023:   DMDAVecGetArray(stokes_da,X_local,&stokes);

1025:   DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1026:   DMDAGetBoundingBox(stokes_da,xymin,xymax);

1028:   h = (xymax[0]-xymin[0])/((PetscReal)M);

1030:   tp_L2 = tu_L2 = tu_H1 = 0.0;

1032:   DMDAGetElementCorners(stokes_da,&sex,&sey,NULL,&mx,&my,NULL);
1033:   for (ej = sey; ej < sey+my; ej++) {
1034:     for (ei = sex; ei < sex+mx; ei++) {
1035:       /* get coords for the element */
1036:       GetElementCoords(_coords,ei,ej,el_coords);
1037:       StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1038:       StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);

1040:       /* evaluate integral */
1041:       p_e_L2 = 0.0;
1042:       u_e_L2 = 0.0;
1043:       u_e_H1 = 0.0;
1044:       for (p = 0; p < ngp; p++) {
1045:         ConstructQ12D_Ni(gp_xi[p],Ni_p);
1046:         ConstructQ12D_GNi(gp_xi[p],GNi_p);
1047:         ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1048:         fac = gp_weight[p]*J_p;

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

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

1055:           u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1056:           v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1057:           u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);

1059:           u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error              /* du/dx */
1060:                                +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error               /* du/dy */
1061:                                +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error               /* dv/dx */
1062:                                +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error);             /* dv/dy */
1063:         }
1064:       }

1066:       tp_L2 += p_e_L2;
1067:       tu_L2 += u_e_L2;
1068:       tu_H1 += u_e_H1;
1069:     }
1070:   }
1071:   MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1072:   MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1073:   MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1074:   p_L2 = PetscSqrtScalar(p_L2);
1075:   u_L2 = PetscSqrtScalar(u_L2);
1076:   u_H1 = PetscSqrtScalar(u_H1);

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

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

1082:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1083:   VecDestroy(&X_analytic_local);
1084:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1085:   VecDestroy(&X_local);
1086:   return(0);
1087: }

1089: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1090: {
1091:   DM                     da_Stokes,da_prop;
1092:   PetscInt               u_dof,p_dof,dof,stencil_width;
1093:   Mat                    A,B;
1094:   PetscInt               mxl,myl;
1095:   DM                     prop_cda,vel_cda;
1096:   Vec                    prop_coords,vel_coords;
1097:   PetscInt               si,sj,nx,ny,i,j,p;
1098:   Vec                    f,X;
1099:   PetscInt               prop_dof,prop_stencil_width;
1100:   Vec                    properties,l_properties;
1101:   PetscReal              dx,dy;
1102:   PetscInt               M,N;
1103:   DMDACoor2d             **_prop_coords,**_vel_coords;
1104:   GaussPointCoefficients **element_props;
1105:   PetscInt               its;
1106:   KSP                    ksp_S;
1107:   PetscInt               coefficient_structure = 0;
1108:   PetscInt               cpu_x,cpu_y,*lx = NULL,*ly = NULL;
1109:   PetscBool              use_gp_coords = PETSC_FALSE,set,output_gnuplot = PETSC_FALSE;
1110:   char                   filename[PETSC_MAX_PATH_LEN];
1111:   PetscErrorCode         ierr;


1115:   PetscOptionsGetBool(NULL,NULL,"-gnuplot",&output_gnuplot,NULL);
1116: 
1117:   /* Generate the da for velocity and pressure */
1118:   /*
1119:   We use Q1 elements for the temperature.
1120:   FEM has a 9-point stencil (BOX) or connectivity pattern
1121:   Num nodes in each direction is mx+1, my+1
1122:   */
1123:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1124:   p_dof         = P_DOFS; /* p - pressure */
1125:   dof           = u_dof+p_dof;
1126:   stencil_width = 1;
1127:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&da_Stokes);
1128:   DMSetFromOptions(da_Stokes);
1129:   DMSetUp(da_Stokes);
1130:   DMDASetFieldName(da_Stokes,0,"Vx");
1131:   DMDASetFieldName(da_Stokes,1,"Vy");
1132:   DMDASetFieldName(da_Stokes,2,"P");

1134:   /* unit box [0,1] x [0,1] */
1135:   DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.,0.);

1137:   /* Generate element properties, we will assume all material properties are constant over the element */
1138:   /* local number of elements */
1139:   DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,NULL);

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

1145:   prop_dof           = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1146:   prop_stencil_width = 0;
1147:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
1148:   DMSetFromOptions(da_prop);
1149:   DMSetUp(da_prop);
1150:   PetscFree(lx);
1151:   PetscFree(ly);

1153:   /* define centroid positions */
1154:   DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1155:   dx   = 1.0/((PetscReal)(M));
1156:   dy   = 1.0/((PetscReal)(N));

1158:   DMDASetUniformCoordinates(da_prop,0.0+0.5*dx,1.0-0.5*dx,0.0+0.5*dy,1.0-0.5*dy,0.,0);

1160:   /* define coefficients */
1161:   PetscOptionsGetInt(NULL,NULL,"-c_str",&coefficient_structure,NULL);

1163:   DMCreateGlobalVector(da_prop,&properties);
1164:   DMCreateLocalVector(da_prop,&l_properties);
1165:   DMDAVecGetArray(da_prop,l_properties,&element_props);

1167:   DMGetCoordinateDM(da_prop,&prop_cda);
1168:   DMGetCoordinatesLocal(da_prop,&prop_coords);
1169:   DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);

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

1173:   DMGetCoordinateDM(da_Stokes,&vel_cda);
1174:   DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1175:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);

1177:   /* interpolate the coordinates */
1178:   for (j = sj; j < sj+ny; j++) {
1179:     for (i = si; i < si+nx; i++) {
1180:       PetscInt    ngp;
1181:       PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1182:       PetscScalar el_coords[8];

1184:       GetElementCoords(_vel_coords,i,j,el_coords);
1185:       ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

1187:       for (p = 0; p < GAUSS_POINTS; p++) {
1188:         PetscScalar gp_x,gp_y;
1189:         PetscInt    n;
1190:         PetscScalar xi_p[2],Ni_p[4];

1192:         xi_p[0] = gp_xi[p][0];
1193:         xi_p[1] = gp_xi[p][1];
1194:         ConstructQ12D_Ni(xi_p,Ni_p);

1196:         gp_x = 0.0;
1197:         gp_y = 0.0;
1198:         for (n = 0; n < NODES_PER_EL; n++) {
1199:           gp_x = gp_x+Ni_p[n]*el_coords[2*n];
1200:           gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1201:         }
1202:         element_props[j][i].gp_coords[2*p]   = gp_x;
1203:         element_props[j][i].gp_coords[2*p+1] = gp_y;
1204:       }
1205:     }
1206:   }

1208:   /* define the coefficients */
1209:   PetscOptionsGetBool(NULL,NULL,"-use_gp_coords",&use_gp_coords,NULL);

1211:   for (j = sj; j < sj+ny; j++) {
1212:     for (i = si; i < si+nx; i++) {
1213:       PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1214:       PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1215:       PetscReal coord_x,coord_y;

1217:       if (coefficient_structure == 0) {
1218:         PetscReal opts_eta0,opts_eta1,opts_xc;
1219:         PetscInt  opts_nz;

1221:         opts_eta0 = 1.0;
1222:         opts_eta1 = 1.0;
1223:         opts_xc   = 0.5;
1224:         opts_nz   = 1;

1226:         PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1227:         PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1228:         PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1229:         PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);

1231:         for (p = 0; p < GAUSS_POINTS; p++) {
1232:           coord_x = centroid_x;
1233:           coord_y = centroid_y;
1234:           if (use_gp_coords) {
1235:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1236:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1237:           }

1239:           element_props[j][i].eta[p] = opts_eta0;
1240:           if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;

1242:           element_props[j][i].fx[p] = 0.0;
1243:           element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1244:         }
1245:       } else if (coefficient_structure == 1) { /* square sinker */
1246:         PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;

1248:         opts_eta0 = 1.0;
1249:         opts_eta1 = 1.0;
1250:         opts_dx   = 0.50;
1251:         opts_dy   = 0.50;

1253:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1254:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1255:         PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1256:         PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);

1258:         for (p = 0; p < GAUSS_POINTS; p++) {
1259:           coord_x = centroid_x;
1260:           coord_y = centroid_y;
1261:           if (use_gp_coords) {
1262:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1263:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1264:           }

1266:           element_props[j][i].eta[p] = opts_eta0;
1267:           element_props[j][i].fx[p]  = 0.0;
1268:           element_props[j][i].fy[p]  = 0.0;

1270:           if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1271:             if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1272:               element_props[j][i].eta[p] =  opts_eta1;
1273:               element_props[j][i].fx[p]  =  0.0;
1274:               element_props[j][i].fy[p]  = -1.0;
1275:             }
1276:           }
1277:         }
1278:       } else if (coefficient_structure == 2) { /* circular sinker */
1279:         PetscReal opts_eta0,opts_eta1,opts_r,radius2;

1281:         opts_eta0 = 1.0;
1282:         opts_eta1 = 1.0;
1283:         opts_r    = 0.25;

1285:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1286:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1287:         PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);

1289:         for (p = 0; p < GAUSS_POINTS; p++) {
1290:           coord_x = centroid_x;
1291:           coord_y = centroid_y;
1292:           if (use_gp_coords) {
1293:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1294:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1295:           }

1297:           element_props[j][i].eta[p] = opts_eta0;
1298:           element_props[j][i].fx[p]  = 0.0;
1299:           element_props[j][i].fy[p]  = 0.0;

1301:           radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1302:           if (radius2 < opts_r*opts_r) {
1303:             element_props[j][i].eta[p] =  opts_eta1;
1304:             element_props[j][i].fx[p]  =  0.0;
1305:             element_props[j][i].fy[p]  = -1.0;
1306:           }
1307:         }
1308:       } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1309:         PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;

1311:         opts_eta0 = 1.0;
1312:         opts_eta1 = 1.0;
1313:         opts_r    = 0.25;
1314:         opts_c0x  = 0.35;       /* circle center */
1315:         opts_c0y  = 0.35;
1316:         opts_s0x  = 0.7;       /* square center */
1317:         opts_s0y  = 0.7;
1318:         opts_dx   = 0.25;
1319:         opts_dy   = 0.25;
1320:         opts_phi  = 25;

1322:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1323:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1324:         PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);
1325:         PetscOptionsGetReal(NULL,NULL,"-sinker_c0x",&opts_c0x,NULL);
1326:         PetscOptionsGetReal(NULL,NULL,"-sinker_c0y",&opts_c0y,NULL);
1327:         PetscOptionsGetReal(NULL,NULL,"-sinker_s0x",&opts_s0x,NULL);
1328:         PetscOptionsGetReal(NULL,NULL,"-sinker_s0y",&opts_s0y,NULL);
1329:         PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1330:         PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
1331:         PetscOptionsGetReal(NULL,NULL,"-sinker_phi",&opts_phi,NULL);
1332:         opts_phi *= PETSC_PI / 180;

1334:         for (p = 0; p < GAUSS_POINTS; p++) {
1335:           coord_x = centroid_x;
1336:           coord_y = centroid_y;
1337:           if (use_gp_coords) {
1338:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1339:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1340:           }

1342:           element_props[j][i].eta[p] = opts_eta0;
1343:           element_props[j][i].fx[p]  = 0.0;
1344:           element_props[j][i].fy[p]  = -0.2;

1346:           radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1347:           if (radius2 < opts_r*opts_r
1348:               || (PetscAbs(+(coord_x - opts_s0x)*PetscCosReal(opts_phi) + (coord_y - opts_s0y)*PetscSinReal(opts_phi)) < opts_dx/2
1349:                   && PetscAbs(-(coord_x - opts_s0x)*PetscSinReal(opts_phi) + (coord_y - opts_s0y)*PetscCosReal(opts_phi)) < opts_dy/2)) {
1350:             element_props[j][i].eta[p] =  opts_eta1;
1351:             element_props[j][i].fx[p]  =  0.0;
1352:             element_props[j][i].fy[p]  = -1.0;
1353:           }
1354:         }
1355:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1356:     }
1357:   }
1358:   DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);

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

1362:   DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1363:   DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1364:   DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);

1366:   if (output_gnuplot) {
1367:     DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1368:     DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for Stokes eqn.","properties");
1369:   }

1371:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1372:   DMSetMatType(da_Stokes,MATAIJ);
1373:   DMCreateMatrix(da_Stokes,&A);
1374:   DMCreateMatrix(da_Stokes,&B);
1375:   DMCreateGlobalVector(da_Stokes,&f);
1376:   DMCreateGlobalVector(da_Stokes,&X);

1378:   /* assemble A11 */
1379:   MatZeroEntries(A);
1380:   MatZeroEntries(B);
1381:   VecZeroEntries(f);

1383:   AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1384:   AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1385:   /* build force vector */
1386:   AssembleF_Stokes(f,da_Stokes,da_prop,properties);

1388:   DMDABCApplyFreeSlip(da_Stokes,A,f);
1389:   DMDABCApplyFreeSlip(da_Stokes,B,NULL);

1391:   /* SOLVE */
1392:   KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1393:   KSPSetOptionsPrefix(ksp_S,"stokes_");
1394:   KSPSetDM(ksp_S,da_Stokes);
1395:   KSPSetDMActive(ksp_S,PETSC_FALSE);
1396:   KSPSetOperators(ksp_S,A,B);
1397:   KSPSetFromOptions(ksp_S);
1398:   {
1399:     PC             pc;
1400:     const PetscInt ufields[] = {0,1},pfields[1] = {2};
1401:     KSPGetPC(ksp_S,&pc);
1402:     PCFieldSplitSetBlockSize(pc,3);
1403:     PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1404:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1405:   }

1407:   {
1408:     PC        pc_S;
1409:     KSP       *sub_ksp,ksp_U;
1410:     PetscInt  nsplits;
1411:     DM        da_U;
1412:     PetscBool is_pcfs;

1414:     KSPSetUp(ksp_S);
1415:     KSPGetPC(ksp_S,&pc_S);

1417:     is_pcfs = PETSC_FALSE;
1418:     PetscObjectTypeCompare((PetscObject)pc_S,PCFIELDSPLIT,&is_pcfs);

1420:     if (is_pcfs) {
1421:       PCFieldSplitGetSubKSP(pc_S,&nsplits,&sub_ksp);
1422:       ksp_U = sub_ksp[0];
1423:       PetscFree(sub_ksp);

1425:       if (nsplits == 2) {
1426:         DMDAGetReducedDMDA(da_Stokes,2,&da_U);

1428:         KSPSetDM(ksp_U,da_U);
1429:         KSPSetDMActive(ksp_U,PETSC_FALSE);
1430:         DMDestroy(&da_U);
1431:       }
1432:     }
1433:   }

1435:   KSPSolve(ksp_S,f,X);

1437:   PetscOptionsGetString(NULL,NULL,"-o",filename,sizeof(filename),&set);
1438:   if (set) {
1439:     char        *ext;
1440:     PetscViewer viewer;
1441:     PetscBool   flg;
1442:     PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1443:     PetscStrrchr(filename,'.',&ext);
1444:     PetscStrcmp("vts",ext,&flg);
1445:     if (flg) {
1446:       PetscViewerSetType(viewer,PETSCVIEWERVTK);
1447:     } else {
1448:       PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1449:     }
1450:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1451:     PetscViewerFileSetName(viewer,filename);
1452:     VecView(X,viewer);
1453:     PetscViewerDestroy(&viewer);
1454:   }
1455:   if (output_gnuplot) {
1456:     DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
1457:   }

1459:   KSPGetIterationNumber(ksp_S,&its);

1461:   if (coefficient_structure == 0) {
1462:     PetscReal opts_eta0,opts_eta1,opts_xc;
1463:     PetscInt  opts_nz,N;
1464:     DM        da_Stokes_analytic;
1465:     Vec       X_analytic;
1466:     PetscReal nrm1[3],nrm2[3],nrmI[3];

1468:     opts_eta0 = 1.0;
1469:     opts_eta1 = 1.0;
1470:     opts_xc   = 0.5;
1471:     opts_nz   = 1;

1473:     PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1474:     PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1475:     PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1476:     PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);

1478:     DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1479:     if (output_gnuplot) {
1480:       DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
1481:     }
1482:     DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);

1484:     VecAXPY(X_analytic,-1.0,X);
1485:     VecGetSize(X_analytic,&N);
1486:     N    = N/3;

1488:     VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1489:     VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1490:     VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);

1492:     VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1493:     VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1494:     VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);

1496:     VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1497:     VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1498:     VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);

1500:     DMDestroy(&da_Stokes_analytic);
1501:     VecDestroy(&X_analytic);
1502:   }

1504:   KSPDestroy(&ksp_S);
1505:   VecDestroy(&X);
1506:   VecDestroy(&f);
1507:   MatDestroy(&A);
1508:   MatDestroy(&B);

1510:   DMDestroy(&da_Stokes);
1511:   DMDestroy(&da_prop);

1513:   VecDestroy(&properties);
1514:   VecDestroy(&l_properties);
1515:   return(0);
1516: }

1518: int main(int argc,char **args)
1519: {
1521:   PetscInt       mx,my;

1523:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
1524:   mx   = my = 10;
1525:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1526:   PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1527:   solve_stokes_2d_coupled(mx,my);
1528:   PetscFinalize();
1529:   return ierr;
1530: }

1532: /* -------------------------- helpers for boundary conditions -------------------------------- */
1533: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
1534: {
1535:   DM                     cda;
1536:   Vec                    coords;
1537:   PetscInt               si,sj,nx,ny,i,j;
1538:   PetscInt               M,N;
1539:   DMDACoor2d             **_coords;
1540:   const PetscInt         *g_idx;
1541:   PetscInt               *bc_global_ids;
1542:   PetscScalar            *bc_vals;
1543:   PetscInt               nbcs;
1544:   PetscInt               n_dofs;
1545:   PetscErrorCode         ierr;
1546:   ISLocalToGlobalMapping ltogm;

1549:   DMGetLocalToGlobalMapping(da,&ltogm);
1550:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1552:   DMGetCoordinateDM(da,&cda);
1553:   DMGetCoordinatesLocal(da,&coords);
1554:   DMDAVecGetArray(cda,coords,&_coords);
1555:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1556:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1558:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
1559:   PetscMalloc1(ny*n_dofs,&bc_vals);

1561:   /* init the entries to -1 so VecSetValues will ignore them */
1562:   for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;

1564:   i = nx-1;
1565:   for (j = 0; j < ny; j++) {
1566:     PetscInt local_id;

1568:     local_id = i+j*nx;

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

1572:     bc_vals[j] =  0.0;
1573:   }
1574:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1575:   nbcs = 0;
1576:   if ((si+nx) == (M)) nbcs = ny;

1578:   if (b) {
1579:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1580:     VecAssemblyBegin(b);
1581:     VecAssemblyEnd(b);
1582:   }
1583:   if (A) {
1584:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1585:   }

1587:   PetscFree(bc_vals);
1588:   PetscFree(bc_global_ids);

1590:   DMDAVecRestoreArray(cda,coords,&_coords);
1591:   return(0);
1592: }

1594: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
1595: {
1596:   DM                     cda;
1597:   Vec                    coords;
1598:   PetscInt               si,sj,nx,ny,i,j;
1599:   PetscInt               M,N;
1600:   DMDACoor2d             **_coords;
1601:   const PetscInt         *g_idx;
1602:   PetscInt               *bc_global_ids;
1603:   PetscScalar            *bc_vals;
1604:   PetscInt               nbcs;
1605:   PetscInt               n_dofs;
1606:   PetscErrorCode         ierr;
1607:   ISLocalToGlobalMapping ltogm;

1610:   DMGetLocalToGlobalMapping(da,&ltogm);
1611:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1613:   DMGetCoordinateDM(da,&cda);
1614:   DMGetCoordinatesLocal(da,&coords);
1615:   DMDAVecGetArray(cda,coords,&_coords);
1616:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1617:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1619:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
1620:   PetscMalloc1(ny*n_dofs,&bc_vals);

1622:   /* init the entries to -1 so VecSetValues will ignore them */
1623:   for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;

1625:   i = 0;
1626:   for (j = 0; j < ny; j++) {
1627:     PetscInt local_id;

1629:     local_id = i+j*nx;

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

1633:     bc_vals[j] =  0.0;
1634:   }
1635:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1636:   nbcs = 0;
1637:   if (si == 0) nbcs = ny;

1639:   if (b) {
1640:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1641:     VecAssemblyBegin(b);
1642:     VecAssemblyEnd(b);
1643:   }

1645:   if (A) {
1646:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1647:   }

1649:   PetscFree(bc_vals);
1650:   PetscFree(bc_global_ids);

1652:   DMDAVecRestoreArray(cda,coords,&_coords);
1653:   return(0);
1654: }

1656: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
1657: {
1658:   DM                     cda;
1659:   Vec                    coords;
1660:   PetscInt               si,sj,nx,ny,i,j;
1661:   PetscInt               M,N;
1662:   DMDACoor2d             **_coords;
1663:   const PetscInt         *g_idx;
1664:   PetscInt               *bc_global_ids;
1665:   PetscScalar            *bc_vals;
1666:   PetscInt               nbcs;
1667:   PetscInt               n_dofs;
1668:   PetscErrorCode         ierr;
1669:   ISLocalToGlobalMapping ltogm;

1672:   DMGetLocalToGlobalMapping(da,&ltogm);
1673:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1675:   DMGetCoordinateDM(da,&cda);
1676:   DMGetCoordinatesLocal(da,&coords);
1677:   DMDAVecGetArray(cda,coords,&_coords);
1678:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1679:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1681:   PetscMalloc1(nx,&bc_global_ids);
1682:   PetscMalloc1(nx,&bc_vals);

1684:   /* init the entries to -1 so VecSetValues will ignore them */
1685:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

1687:   j = ny-1;
1688:   for (i = 0; i < nx; i++) {
1689:     PetscInt local_id;

1691:     local_id = i+j*nx;

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

1695:     bc_vals[i] =  0.0;
1696:   }
1697:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1698:   nbcs = 0;
1699:   if ((sj+ny) == (N)) nbcs = nx;

1701:   if (b) {
1702:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1703:     VecAssemblyBegin(b);
1704:     VecAssemblyEnd(b);
1705:   }
1706:   if (A) {
1707:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1708:   }

1710:   PetscFree(bc_vals);
1711:   PetscFree(bc_global_ids);

1713:   DMDAVecRestoreArray(cda,coords,&_coords);
1714:   return(0);
1715: }

1717: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
1718: {
1719:   DM                     cda;
1720:   Vec                    coords;
1721:   PetscInt               si,sj,nx,ny,i,j;
1722:   PetscInt               M,N;
1723:   DMDACoor2d             **_coords;
1724:   const PetscInt         *g_idx;
1725:   PetscInt               *bc_global_ids;
1726:   PetscScalar            *bc_vals;
1727:   PetscInt               nbcs;
1728:   PetscInt               n_dofs;
1729:   PetscErrorCode         ierr;
1730:   ISLocalToGlobalMapping ltogm;

1733:   DMGetLocalToGlobalMapping(da,&ltogm);
1734:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1736:   DMGetCoordinateDM(da,&cda);
1737:   DMGetCoordinatesLocal(da,&coords);
1738:   DMDAVecGetArray(cda,coords,&_coords);
1739:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1740:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1742:   PetscMalloc1(nx,&bc_global_ids);
1743:   PetscMalloc1(nx,&bc_vals);

1745:   /* init the entries to -1 so VecSetValues will ignore them */
1746:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

1748:   j = 0;
1749:   for (i = 0; i < nx; i++) {
1750:     PetscInt local_id;

1752:     local_id = i+j*nx;

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

1756:     bc_vals[i] =  0.0;
1757:   }
1758:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1759:   nbcs = 0;
1760:   if (sj == 0) nbcs = nx;

1762:   if (b) {
1763:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1764:     VecAssemblyBegin(b);
1765:     VecAssemblyEnd(b);
1766:   }
1767:   if (A) {
1768:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1769:   }

1771:   PetscFree(bc_vals);
1772:   PetscFree(bc_global_ids);

1774:   DMDAVecRestoreArray(cda,coords,&_coords);
1775:   return(0);
1776: }

1778: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
1779: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1780: {

1784:   BCApplyZero_NORTH(da_Stokes,1,A,f);
1785:   BCApplyZero_EAST(da_Stokes,0,A,f);
1786:   BCApplyZero_SOUTH(da_Stokes,1,A,f);
1787:   BCApplyZero_WEST(da_Stokes,0,A,f);
1788:   return(0);
1789: }