Actual source code: ex43.c

  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:      -glvis : Visualizes coefficients and u,p solution through GLVIs (use -viewer_glvis_dmda_bs 2,1 to visualize velocity as a vector)\n\
 11:      -c_str : Indicates the structure of the coefficients to use \n"
 12: "\
 13:           -c_str 0 => Coefficient definition for an analytic solution with a vertical jump in viscosity at x = xc \n\
 14:                       This problem is driven by the forcing function f(x,y) = (0, sin(nz pi y)cos(pi x) \n\
 15:                          Parameters: \n\
 16:                               -solcx_eta0  : Viscosity to the left of the interface \n\
 17:                               -solcx_eta1  : Viscosity to the right of the interface \n\
 18:                               -solcx_xc    : Location of the interface \n\
 19:                               -solcx_nz    : Wavenumber in the y direction \n"
 20: "\
 21:           -c_str 1 => Coefficient definition for a dense rectangular blob located at the center of the domain \n\
 22:                          Parameters: \n\
 23:                               -sinker_eta0 : Viscosity of the background fluid \n\
 24:                               -sinker_eta1 : Viscosity of the blob \n\
 25:                               -sinker_dx   : Width of the blob \n\
 26:                               -sinker_dy   : Height of the blob \n"
 27: "\
 28:           -c_str 2 => Coefficient definition for a dense circular blob located at the center of the domain \n\
 29:                          Parameters: \n\
 30:                               -sinker_eta0 : Viscosity of the background fluid \n\
 31:                               -sinker_eta1 : Viscosity of the blob \n\
 32:                               -sinker_r    : Radius of the blob \n"
 33: "\
 34:           -c_str 3 => Coefficient definition for a dense circular and rectangular inclusion (located at the center of the domain) \n\
 35:                               -sinker_eta0 : Viscosity of the background fluid \n\
 36:                               -sinker_eta1 : Viscosity of the two inclusions \n\
 37:                               -sinker_r    : Radius of the circular inclusion \n\
 38:                               -sinker_c0x  : Origin (x-coord) of the circular inclusion \n\
 39:                               -sinker_c0y  : Origin (y-coord) of the circular inclusion \n\
 40:                               -sinker_dx   : Width of the rectangular inclusion \n\
 41:                               -sinker_dy   : Height of the rectangular inclusion \n\
 42:                               -sinker_phi  : Rotation angle of the rectangular inclusion \n"
 43: "\
 44:           -c_str 4 => Coefficient definition for checkerboard jumps aligned with the domain decomposition \n\
 45:                               -jump_eta0      : Viscosity for black subdomains \n\
 46:                               -jump_magnitude : Magnitude of jumps. White subdomains will have eta = eta0*10^magnitude \n\
 47:                               -jump_nz        : Wavenumber in the y direction for rhs \n"
 48: "\
 49:      -use_gp_coords : Evaluate the viscosity and force term at the global coordinates of each quadrature point \n\
 50:                       By default, the viscosity and force term are evaulated at the element center and applied as a constant over the entire element \n";

 52: /* Contributed by Dave May */

 54: #include <petscksp.h>
 55: #include <petscdm.h>
 56: #include <petscdmda.h>

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

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

 63: #define NSD            2 /* number of spatial dimensions */
 64: #define NODES_PER_EL   4 /* nodes per element */
 65: #define U_DOFS         2 /* degrees of freedom per velocity node */
 66: #define P_DOFS         1 /* degrees of freedom per pressure node */
 67: #define GAUSS_POINTS   4

 69: /* Gauss point based evaluation 8+4+4+4 = 20 */
 70: typedef struct {
 71:   PetscScalar gp_coords[2*GAUSS_POINTS];
 72:   PetscScalar eta[GAUSS_POINTS];
 73:   PetscScalar fx[GAUSS_POINTS];
 74:   PetscScalar fy[GAUSS_POINTS];
 75: } GaussPointCoefficients;

 77: typedef struct {
 78:   PetscScalar u_dof;
 79:   PetscScalar v_dof;
 80:   PetscScalar p_dof;
 81: } StokesDOF;

 83: static PetscErrorCode glvis_extract_eta(PetscObject oV,PetscInt nf, PetscObject oVf[], void *ctx)
 84: {
 85:   DM                     properties_da = (DM)(ctx),stokes_da;
 86:   Vec                    V = (Vec)oV, *Vf = (Vec*)oVf;
 87:   GaussPointCoefficients **props;
 88:   PetscInt               sex,sey,mx,my;
 89:   PetscInt               ei,ej,p,cum;
 90:   PetscScalar            *array;
 91:   PetscErrorCode         ierr;

 94:   VecGetDM(Vf[0],&stokes_da);
 95:   DMDAVecGetArrayRead(properties_da,V,&props);
 96:   DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
 97:   DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
 98:   VecGetArray(Vf[0],&array);
 99:   cum  = 0;
100:   for (ej = sey; ej < sey+my; ej++) {
101:     for (ei = sex; ei < sex+mx; ei++) {
102:       for (p = 0; p < GAUSS_POINTS; p++) {
103:         array[cum++] = props[ej][ei].eta[p];
104:       }
105:     }
106:   }
107:   VecRestoreArray(Vf[0],&array);
108:   DMDAVecRestoreArrayRead(properties_da,V,&props);
109:   return(0);
110: }

112: /*
113:  Element: Local basis function ordering
114:  1-----2
115:  |     |
116:  |     |
117:  0-----3
118:  */
119: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
120: {
121:   PetscScalar xi  = _xi[0];
122:   PetscScalar eta = _xi[1];

124:   Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
125:   Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
126:   Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
127:   Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
128: }

130: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
131: {
132:   PetscScalar xi  = _xi[0];
133:   PetscScalar eta = _xi[1];

135:   GNi[0][0] = -0.25*(1.0-eta);
136:   GNi[0][1] = -0.25*(1.0+eta);
137:   GNi[0][2] =   0.25*(1.0+eta);
138:   GNi[0][3] =   0.25*(1.0-eta);

140:   GNi[1][0] = -0.25*(1.0-xi);
141:   GNi[1][1] =   0.25*(1.0-xi);
142:   GNi[1][2] =   0.25*(1.0+xi);
143:   GNi[1][3] = -0.25*(1.0+xi);
144: }

146: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
147: {
148:   PetscScalar J00,J01,J10,J11,J;
149:   PetscScalar iJ00,iJ01,iJ10,iJ11;
150:   PetscInt    i;

152:   J00 = J01 = J10 = J11 = 0.0;
153:   for (i = 0; i < NODES_PER_EL; i++) {
154:     PetscScalar cx = coords[2*i];
155:     PetscScalar cy = coords[2*i+1];

157:     J00 = J00+GNi[0][i]*cx;      /* J_xx = dx/dxi */
158:     J01 = J01+GNi[0][i]*cy;      /* J_xy = dy/dxi */
159:     J10 = J10+GNi[1][i]*cx;      /* J_yx = dx/deta */
160:     J11 = J11+GNi[1][i]*cy;      /* J_yy = dy/deta */
161:   }
162:   J = (J00*J11)-(J01*J10);

164:   iJ00 =  J11/J;
165:   iJ01 = -J01/J;
166:   iJ10 = -J10/J;
167:   iJ11 =  J00/J;

169:   for (i = 0; i < NODES_PER_EL; i++) {
170:     GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
171:     GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
172:   }

174:   if (det_J) *det_J = J;
175: }

177: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
178: {
179:   *ngp         = 4;
180:   gp_xi[0][0]  = -0.57735026919; gp_xi[0][1] = -0.57735026919;
181:   gp_xi[1][0]  = -0.57735026919; gp_xi[1][1] =  0.57735026919;
182:   gp_xi[2][0]  =  0.57735026919; gp_xi[2][1] =  0.57735026919;
183:   gp_xi[3][0]  =  0.57735026919; gp_xi[3][1] = -0.57735026919;
184:   gp_weight[0] = 1.0;
185:   gp_weight[1] = 1.0;
186:   gp_weight[2] = 1.0;
187:   gp_weight[3] = 1.0;
188: }

190: /*
191: i,j are the element indices
192: The unknown is a vector quantity.
193: The s[].c is used to indicate the degree of freedom.
194: */
195: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
196: {
198:   /* velocity */
199:   /* node 0 */
200:   s_u[0].i = i; s_u[0].j = j; s_u[0].c = 0;       /* Vx0 */
201:   s_u[1].i = i; s_u[1].j = j; s_u[1].c = 1;       /* Vy0 */

203:   /* node 1 */
204:   s_u[2].i = i; s_u[2].j = j+1; s_u[2].c = 0;     /* Vx1 */
205:   s_u[3].i = i; s_u[3].j = j+1; s_u[3].c = 1;     /* Vy1 */

207:   /* node 2 */
208:   s_u[4].i = i+1; s_u[4].j = j+1; s_u[4].c = 0;   /* Vx2 */
209:   s_u[5].i = i+1; s_u[5].j = j+1; s_u[5].c = 1;   /* Vy2 */

211:   /* node 3 */
212:   s_u[6].i = i+1; s_u[6].j = j; s_u[6].c = 0;     /* Vx3 */
213:   s_u[7].i = i+1; s_u[7].j = j; s_u[7].c = 1;     /* Vy3 */

215:   /* pressure */
216:   s_p[0].i = i;   s_p[0].j = j;   s_p[0].c = 2; /* P0 */
217:   s_p[1].i = i;   s_p[1].j = j+1; s_p[1].c = 2; /* P1 */
218:   s_p[2].i = i+1; s_p[2].j = j+1; s_p[2].c = 2; /* P2 */
219:   s_p[3].i = i+1; s_p[3].j = j;   s_p[3].c = 2; /* P3 */
220:   return(0);
221: }

223: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
224: {
226:   PetscMPIInt    rank;
227:   PetscInt       proc_I,proc_J;
228:   PetscInt       cpu_x,cpu_y;
229:   PetscInt       local_mx,local_my;
230:   Vec            vlx,vly;
231:   PetscInt       *LX,*LY,i;
232:   PetscScalar    *_a;
233:   Vec            V_SEQ;
234:   VecScatter     ctx;

237:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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

241:   proc_J = rank/cpu_x;
242:   proc_I = rank-cpu_x*proc_J;

244:   PetscMalloc1(cpu_x,&LX);
245:   PetscMalloc1(cpu_y,&LY);

247:   DMDAGetElementsSizes(da,&local_mx,&local_my,NULL);
248:   VecCreate(PETSC_COMM_WORLD,&vlx);
249:   VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
250:   VecSetFromOptions(vlx);

252:   VecCreate(PETSC_COMM_WORLD,&vly);
253:   VecSetSizes(vly,PETSC_DECIDE,cpu_y);
254:   VecSetFromOptions(vly);

256:   VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
257:   VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
258:   VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
259:   VecAssemblyBegin(vly);VecAssemblyEnd(vly);

261:   VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
262:   VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
263:   VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
264:   VecGetArray(V_SEQ,&_a);
265:   for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
266:   VecRestoreArray(V_SEQ,&_a);
267:   VecScatterDestroy(&ctx);
268:   VecDestroy(&V_SEQ);

270:   VecScatterCreateToAll(vly,&ctx,&V_SEQ);
271:   VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
272:   VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
273:   VecGetArray(V_SEQ,&_a);
274:   for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
275:   VecRestoreArray(V_SEQ,&_a);
276:   VecScatterDestroy(&ctx);
277:   VecDestroy(&V_SEQ);

279:   *_lx = LX;
280:   *_ly = LY;

282:   VecDestroy(&vlx);
283:   VecDestroy(&vly);
284:   return(0);
285: }

287: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
288: {
289:   DM             cda;
290:   Vec            coords;
291:   DMDACoor2d     **_coords;
292:   PetscInt       si,sj,nx,ny,i,j;
293:   FILE           *fp;
294:   char           fname[PETSC_MAX_PATH_LEN];
295:   PetscMPIInt    rank;

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

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

306:   DMGetCoordinateDM(da,&cda);
307:   DMGetCoordinatesLocal(da,&coords);
308:   DMDAVecGetArrayRead(cda,coords,&_coords);
309:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
310:   for (j = sj; j < sj+ny-1; j++) {
311:     for (i = si; i < si+nx-1; i++) {
312:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
313:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i].x),(double)PetscRealPart(_coords[j+1][i].y));
314:       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));
315:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i+1].x),(double)PetscRealPart(_coords[j][i+1].y));
316:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
317:     }
318:   }
319:   DMDAVecRestoreArrayRead(cda,coords,&_coords);

321:   PetscFClose(PETSC_COMM_SELF,fp);
322:   return(0);
323: }

325: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
326: {
327:   DM             cda;
328:   Vec            coords,local_fields;
329:   DMDACoor2d     **_coords;
330:   FILE           *fp;
331:   char           fname[PETSC_MAX_PATH_LEN];
332:   PetscMPIInt    rank;
333:   PetscInt       si,sj,nx,ny,i,j;
334:   PetscInt       n_dofs,d;
335:   PetscScalar    *_fields;

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

344:   PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
345:   DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
346:   PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
347:   for (d = 0; d < n_dofs; d++) {
348:     const char *field_name;
349:     DMDAGetFieldName(da,d,&field_name);
350:     PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
351:   }
352:   PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");

354:   DMGetCoordinateDM(da,&cda);
355:   DMGetCoordinatesLocal(da,&coords);
356:   DMDAVecGetArray(cda,coords,&_coords);
357:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

359:   DMCreateLocalVector(da,&local_fields);
360:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
361:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
362:   VecGetArray(local_fields,&_fields);

364:   for (j = sj; j < sj+ny; j++) {
365:     for (i = si; i < si+nx; i++) {
366:       PetscScalar coord_x,coord_y;
367:       PetscScalar field_d;

369:       coord_x = _coords[j][i].x;
370:       coord_y = _coords[j][i].y;

372:       PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
373:       for (d = 0; d < n_dofs; d++) {
374:         field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
375:         PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
376:       }
377:       PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
378:     }
379:   }
380:   VecRestoreArray(local_fields,&_fields);
381:   VecDestroy(&local_fields);

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

385:   PetscFClose(PETSC_COMM_SELF,fp);
386:   return(0);
387: }

389: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
390: {
391:   DM                     cda;
392:   Vec                    local_fields;
393:   FILE                   *fp;
394:   char                   fname[PETSC_MAX_PATH_LEN];
395:   PetscMPIInt            rank;
396:   PetscInt               si,sj,nx,ny,i,j,p;
397:   PetscInt               n_dofs,d;
398:   GaussPointCoefficients **_coefficients;
399:   PetscErrorCode         ierr;

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

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

417:   DMGetCoordinateDM(da,&cda);
418:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

420:   DMCreateLocalVector(da,&local_fields);
421:   DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
422:   DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
423:   DMDAVecGetArray(da,local_fields,&_coefficients);

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

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

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

435:         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]));
436:         PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
437:       }
438:     }
439:   }
440:   DMDAVecRestoreArray(da,local_fields,&_coefficients);
441:   VecDestroy(&local_fields);

443:   PetscFClose(PETSC_COMM_SELF,fp);
444:   return(0);
445: }

447: 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)
448: {
449:   PetscInt ij;
450:   PetscInt r,c,nc;

452:   nc = u_NPE*u_dof;
453:   r = w_dof*wi+wd;
454:   c = u_dof*ui+ud;
455:   ij = r*nc+c;
456:   return ij;
457: }

459: /*
460:  D = [ 2.eta   0   0   ]
461:      [   0   2.eta 0   ]
462:      [   0     0   eta ]

464:  B = [ d_dx   0   ]
465:      [  0    d_dy ]
466:      [ d_dy  d_dx ]
467: */
468: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
469: {
470:   PetscInt       ngp;
471:   PetscScalar    gp_xi[GAUSS_POINTS][2];
472:   PetscScalar    gp_weight[GAUSS_POINTS];
473:   PetscInt       p,i,j,k;
474:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
475:   PetscScalar    J_p,tildeD[3];
476:   PetscScalar    B[3][U_DOFS*NODES_PER_EL];

478:   /* define quadrature rule */
479:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

481:   /* evaluate integral */
482:   for (p = 0; p < ngp; p++) {
483:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
484:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);

486:     for (i = 0; i < NODES_PER_EL; i++) {
487:       PetscScalar d_dx_i = GNx_p[0][i];
488:       PetscScalar d_dy_i = GNx_p[1][i];

490:       B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
491:       B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
492:       B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
493:     }

495:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
496:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
497:     tildeD[2] =       gp_weight[p]*J_p*eta[p];

499:     /* form Bt tildeD B */
500:     /*
501:     Ke_ij = Bt_ik . D_kl . B_lj
502:           = B_ki . D_kl . B_lj
503:           = B_ki . D_kk . B_kj
504:     */
505:     for (i = 0; i < 8; i++) {
506:       for (j = 0; j < 8; j++) {
507:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
508:           Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
509:         }
510:       }
511:     }
512:   }
513: }

515: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
516: {
517:   PetscInt    ngp;
518:   PetscScalar gp_xi[GAUSS_POINTS][2];
519:   PetscScalar gp_weight[GAUSS_POINTS];
520:   PetscInt    p,i,j,di;
521:   PetscScalar Ni_p[NODES_PER_EL];
522:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
523:   PetscScalar J_p,fac;

525:   /* define quadrature rule */
526:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

528:   /* evaluate integral */
529:   for (p = 0; p < ngp; p++) {
530:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
531:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
532:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
533:     fac = gp_weight[p]*J_p;

535:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
536:       for (di = 0; di < NSD; di++) { /* u dofs */
537:         for (j = 0; j < 4; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
538:           PetscInt IJ;
539:           IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);

541:           Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
542:         }
543:       }
544:     }
545:   }
546: }

548: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
549: {
550:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
551:   PetscInt    i,j;
552:   PetscInt    nr_g,nc_g;

554:   PetscMemzero(Ge,sizeof(Ge));
555:   FormGradientOperatorQ1(Ge,coords);

557:   nr_g = U_DOFS*NODES_PER_EL;
558:   nc_g = P_DOFS*NODES_PER_EL;

560:   for (i = 0; i < nr_g; i++) {
561:     for (j = 0; j < nc_g; j++) {
562:       De[nr_g*j+i] = Ge[nc_g*i+j];
563:     }
564:   }
565: }

567: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
568: {
569:   PetscInt    ngp;
570:   PetscScalar gp_xi[GAUSS_POINTS][2];
571:   PetscScalar gp_weight[GAUSS_POINTS];
572:   PetscInt    p,i,j;
573:   PetscScalar Ni_p[NODES_PER_EL];
574:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
575:   PetscScalar J_p,fac,eta_avg;

577:   /* define quadrature rule */
578:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

580:   /* evaluate integral */
581:   for (p = 0; p < ngp; p++) {
582:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
583:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
584:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
585:     fac = gp_weight[p]*J_p;

587:     for (i = 0; i < NODES_PER_EL; i++) {
588:       for (j = 0; j < NODES_PER_EL; j++) {
589:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
590:       }
591:     }
592:   }

594:   /* scale */
595:   eta_avg = 0.0;
596:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
597:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
598:   fac     = 1.0/eta_avg;
599:   for (i = 0; i < NODES_PER_EL; i++) {
600:     for (j = 0; j < NODES_PER_EL; j++) {
601:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
602:     }
603:   }
604: }

606: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
607: {
608:   PetscInt    ngp;
609:   PetscScalar gp_xi[GAUSS_POINTS][2];
610:   PetscScalar gp_weight[GAUSS_POINTS];
611:   PetscInt    p,i,j;
612:   PetscScalar Ni_p[NODES_PER_EL];
613:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
614:   PetscScalar J_p,fac,eta_avg;

616:   /* define quadrature rule */
617:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

619:   /* evaluate integral */
620:   for (p = 0; p < ngp; p++) {
621:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
622:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
623:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
624:     fac = gp_weight[p]*J_p;

626:     for (i = 0; i < NODES_PER_EL; i++) {
627:       for (j = 0; j < NODES_PER_EL; j++) {
628:         Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
629:       }
630:     }
631:   }

633:   /* scale */
634:   eta_avg = 0.0;
635:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
636:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
637:   fac     = 1.0/eta_avg;
638:   for (i = 0; i < NODES_PER_EL; i++) {
639:     for (j = 0; j < NODES_PER_EL; j++) {
640:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
641:     }
642:   }
643: }

645: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
646: {
647:   PetscInt    ngp;
648:   PetscScalar gp_xi[GAUSS_POINTS][2];
649:   PetscScalar gp_weight[GAUSS_POINTS];
650:   PetscInt    p,i;
651:   PetscScalar Ni_p[NODES_PER_EL];
652:   PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
653:   PetscScalar J_p,fac;

655:   /* define quadrature rule */
656:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

658:   /* evaluate integral */
659:   for (p = 0; p < ngp; p++) {
660:     ConstructQ12D_Ni(gp_xi[p],Ni_p);
661:     ConstructQ12D_GNi(gp_xi[p],GNi_p);
662:     ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
663:     fac = gp_weight[p]*J_p;

665:     for (i = 0; i < NODES_PER_EL; i++) {
666:       Fe[NSD*i]   += fac*Ni_p[i]*fx[p];
667:       Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
668:     }
669:   }
670: }

672: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
673: {
675:   /* get coords for the element */
676:   el_coords[NSD*0] = _coords[ej][ei].x;     el_coords[NSD*0+1] = _coords[ej][ei].y;
677:   el_coords[NSD*1] = _coords[ej+1][ei].x;   el_coords[NSD*1+1] = _coords[ej+1][ei].y;
678:   el_coords[NSD*2] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
679:   el_coords[NSD*3] = _coords[ej][ei+1].x;   el_coords[NSD*3+1] = _coords[ej][ei+1].y;
680:   return(0);
681: }

683: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
684: {
685:   DM                     cda;
686:   Vec                    coords;
687:   DMDACoor2d             **_coords;
688:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
689:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
690:   PetscInt               sex,sey,mx,my;
691:   PetscInt               ei,ej;
692:   PetscScalar            Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
693:   PetscScalar            Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
694:   PetscScalar            De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
695:   PetscScalar            Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
696:   PetscScalar            el_coords[NODES_PER_EL*NSD];
697:   Vec                    local_properties;
698:   GaussPointCoefficients **props;
699:   PetscScalar            *prop_eta;
700:   PetscErrorCode         ierr;

703:   /* setup for coords */
704:   DMGetCoordinateDM(stokes_da,&cda);
705:   DMGetCoordinatesLocal(stokes_da,&coords);
706:   DMDAVecGetArray(cda,coords,&_coords);

708:   /* setup for coefficients */
709:   DMCreateLocalVector(properties_da,&local_properties);
710:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
711:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
712:   DMDAVecGetArray(properties_da,local_properties,&props);

714:   DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
715:   DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
716:   for (ej = sey; ej < sey+my; ej++) {
717:     for (ei = sex; ei < sex+mx; ei++) {
718:       /* get coords for the element */
719:       GetElementCoords(_coords,ei,ej,el_coords);

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

724:       /* initialise element stiffness matrix */
725:       PetscMemzero(Ae,sizeof(Ae));
726:       PetscMemzero(Ge,sizeof(Ge));
727:       PetscMemzero(De,sizeof(De));
728:       PetscMemzero(Ce,sizeof(Ce));

730:       /* form element stiffness matrix */
731:       FormStressOperatorQ1(Ae,el_coords,prop_eta);
732:       FormGradientOperatorQ1(Ge,el_coords);
733:       FormDivergenceOperatorQ1(De,el_coords);
734:       FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);

736:       /* insert element matrix into global matrix */
737:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
738:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
739:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
740:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
741:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
742:     }
743:   }
744:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
745:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

749:   DMDAVecRestoreArray(properties_da,local_properties,&props);
750:   VecDestroy(&local_properties);
751:   return(0);
752: }

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

774:   /* setup for coords */
775:   DMGetCoordinateDM(stokes_da,&cda);
776:   DMGetCoordinatesLocal(stokes_da,&coords);
777:   DMDAVecGetArray(cda,coords,&_coords);

779:   /* setup for coefficients */
780:   DMCreateLocalVector(properties_da,&local_properties);
781:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
782:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
783:   DMDAVecGetArray(properties_da,local_properties,&props);

785:   DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
786:   DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
787:   for (ej = sey; ej < sey+my; ej++) {
788:     for (ei = sex; ei < sex+mx; ei++) {
789:       /* get coords for the element */
790:       GetElementCoords(_coords,ei,ej,el_coords);

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

795:       /* initialise element stiffness matrix */
796:       PetscMemzero(Ae,sizeof(Ae));
797:       PetscMemzero(Ge,sizeof(Ge));
798:       PetscMemzero(De,sizeof(De));
799:       PetscMemzero(Ce,sizeof(Ce));

801:       /* form element stiffness matrix */
802:       FormStressOperatorQ1(Ae,el_coords,prop_eta);
803:       FormGradientOperatorQ1(Ge,el_coords);
804:       FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);

806:       /* insert element matrix into global matrix */
807:       DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
808:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
809:       MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
810:       MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
811:     }
812:   }
813:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
814:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

818:   DMDAVecRestoreArray(properties_da,local_properties,&props);
819:   VecDestroy(&local_properties);
820:   return(0);
821: }

823: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
824: {
825:   PetscInt n;

828:   for (n = 0; n < 4; n++) {
829:     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];
830:     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];
831:     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];
832:   }
833:   return(0);
834: }

836: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
837: {
838:   DM                     cda;
839:   Vec                    coords;
840:   DMDACoor2d             **_coords;
841:   MatStencil             u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
842:   MatStencil             p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
843:   PetscInt               sex,sey,mx,my;
844:   PetscInt               ei,ej;
845:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
846:   PetscScalar            He[NODES_PER_EL*P_DOFS];
847:   PetscScalar            el_coords[NODES_PER_EL*NSD];
848:   Vec                    local_properties;
849:   GaussPointCoefficients **props;
850:   PetscScalar            *prop_fx,*prop_fy;
851:   Vec                    local_F;
852:   StokesDOF              **ff;
853:   PetscErrorCode         ierr;

856:   /* setup for coords */
857:   DMGetCoordinateDM(stokes_da,&cda);
858:   DMGetCoordinatesLocal(stokes_da,&coords);
859:   DMDAVecGetArray(cda,coords,&_coords);

861:   /* setup for coefficients */
862:   DMGetLocalVector(properties_da,&local_properties);
863:   DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
864:   DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
865:   DMDAVecGetArray(properties_da,local_properties,&props);

867:   /* get access to the vector */
868:   DMGetLocalVector(stokes_da,&local_F);
869:   VecZeroEntries(local_F);
870:   DMDAVecGetArray(stokes_da,local_F,&ff);

872:   DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
873:   DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
874:   for (ej = sey; ej < sey+my; ej++) {
875:     for (ei = sex; ei < sex+mx; ei++) {
876:       /* get coords for the element */
877:       GetElementCoords(_coords,ei,ej,el_coords);

879:       /* get coefficients for the element */
880:       prop_fx = props[ej][ei].fx;
881:       prop_fy = props[ej][ei].fy;

883:       /* initialise element stiffness matrix */
884:       PetscMemzero(Fe,sizeof(Fe));
885:       PetscMemzero(He,sizeof(He));

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

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

893:       DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
894:     }
895:   }

897:   DMDAVecRestoreArray(stokes_da,local_F,&ff);
898:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
899:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
900:   DMRestoreLocalVector(stokes_da,&local_F);

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

904:   DMDAVecRestoreArray(properties_da,local_properties,&props);
905:   DMRestoreLocalVector(properties_da,&local_properties);
906:   return(0);
907: }

909: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM *_da,Vec *_X)
910: {
911:   DM             da,cda;
912:   Vec            X;
913:   StokesDOF      **_stokes;
914:   Vec            coords;
915:   DMDACoor2d     **_coords;
916:   PetscInt       si,sj,ei,ej,i,j;

920:   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);
921:   DMSetFromOptions(da);
922:   DMSetUp(da);
923:   DMDASetFieldName(da,0,"anlytic_Vx");
924:   DMDASetFieldName(da,1,"anlytic_Vy");
925:   DMDASetFieldName(da,2,"analytic_P");

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

929:   DMGetCoordinatesLocal(da,&coords);
930:   DMGetCoordinateDM(da,&cda);
931:   DMDAVecGetArray(cda,coords,&_coords);

933:   DMCreateGlobalVector(da,&X);
934:   DMDAVecGetArray(da,X,&_stokes);

936:   DMDAGetCorners(da,&si,&sj,0,&ei,&ej,0);
937:   for (j = sj; j < sj+ej; j++) {
938:     for (i = si; i < si+ei; i++) {
939:       PetscReal pos[2],pressure,vel[2],total_stress[3],strain_rate[3];

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

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

946:       _stokes[j][i].u_dof = vel[0];
947:       _stokes[j][i].v_dof = vel[1];
948:       _stokes[j][i].p_dof = pressure;
949:     }
950:   }
951:   DMDAVecRestoreArray(da,X,&_stokes);
952:   DMDAVecRestoreArray(cda,coords,&_coords);

954:   *_da = da;
955:   *_X  = X;
956:   return(0);
957: }

959: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
960: {
962:   /* get the nodal fields */
963:   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;
964:   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;
965:   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;
966:   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;
967:   return(0);
968: }

970: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
971: {
972:   DM             cda;
973:   Vec            coords,X_analytic_local,X_local;
974:   DMDACoor2d     **_coords;
975:   PetscInt       sex,sey,mx,my;
976:   PetscInt       ei,ej;
977:   PetscScalar    el_coords[NODES_PER_EL*NSD];
978:   StokesDOF      **stokes_analytic,**stokes;
979:   StokesDOF      stokes_analytic_e[4],stokes_e[4];

981:   PetscScalar    GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
982:   PetscScalar    Ni_p[NODES_PER_EL];
983:   PetscInt       ngp;
984:   PetscScalar    gp_xi[GAUSS_POINTS][2];
985:   PetscScalar    gp_weight[GAUSS_POINTS];
986:   PetscInt       p,i;
987:   PetscScalar    J_p,fac;
988:   PetscScalar    h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
989:   PetscInt       M;
990:   PetscReal      xymin[2],xymax[2];

994:   /* define quadrature rule */
995:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

997:   /* setup for coords */
998:   DMGetCoordinateDM(stokes_da,&cda);
999:   DMGetCoordinatesLocal(stokes_da,&coords);
1000:   DMDAVecGetArray(cda,coords,&_coords);

1002:   /* setup for analytic */
1003:   DMCreateLocalVector(stokes_da,&X_analytic_local);
1004:   DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1005:   DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1006:   DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);

1008:   /* setup for solution */
1009:   DMCreateLocalVector(stokes_da,&X_local);
1010:   DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1011:   DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1012:   DMDAVecGetArray(stokes_da,X_local,&stokes);

1014:   DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1015:   DMGetBoundingBox(stokes_da,xymin,xymax);

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

1019:   tp_L2 = tu_L2 = tu_H1 = 0.0;

1021:   DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
1022:   DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
1023:   for (ej = sey; ej < sey+my; ej++) {
1024:     for (ei = sex; ei < sex+mx; ei++) {
1025:       /* get coords for the element */
1026:       GetElementCoords(_coords,ei,ej,el_coords);
1027:       StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1028:       StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);

1030:       /* evaluate integral */
1031:       p_e_L2 = 0.0;
1032:       u_e_L2 = 0.0;
1033:       u_e_H1 = 0.0;
1034:       for (p = 0; p < ngp; p++) {
1035:         ConstructQ12D_Ni(gp_xi[p],Ni_p);
1036:         ConstructQ12D_GNi(gp_xi[p],GNi_p);
1037:         ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1038:         fac = gp_weight[p]*J_p;

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

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

1045:           u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1046:           v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1047:           u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);

1049:           u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error              /* du/dx */
1050:                                +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error               /* du/dy */
1051:                                +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error               /* dv/dx */
1052:                                +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error);             /* dv/dy */
1053:         }
1054:       }

1056:       tp_L2 += p_e_L2;
1057:       tu_L2 += u_e_L2;
1058:       tu_H1 += u_e_H1;
1059:     }
1060:   }
1061:   MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1062:   MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1063:   MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1064:   p_L2 = PetscSqrtScalar(p_L2);
1065:   u_L2 = PetscSqrtScalar(u_L2);
1066:   u_H1 = PetscSqrtScalar(u_H1);

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

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

1072:   DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1073:   VecDestroy(&X_analytic_local);
1074:   DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1075:   VecDestroy(&X_local);
1076:   return(0);
1077: }

1079: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1080: {
1081:   DM                     da_Stokes,da_prop;
1082:   PetscInt               u_dof,p_dof,dof,stencil_width;
1083:   Mat                    A,B;
1084:   DM                     prop_cda,vel_cda;
1085:   Vec                    prop_coords,vel_coords;
1086:   PetscInt               si,sj,nx,ny,i,j,p;
1087:   Vec                    f,X;
1088:   PetscInt               prop_dof,prop_stencil_width;
1089:   Vec                    properties,l_properties;
1090:   PetscReal              dx,dy;
1091:   PetscInt               M,N;
1092:   DMDACoor2d             **_prop_coords,**_vel_coords;
1093:   GaussPointCoefficients **element_props;
1094:   PetscInt               its;
1095:   KSP                    ksp_S;
1096:   PetscInt               coefficient_structure = 0;
1097:   PetscInt               cpu_x,cpu_y,*lx = NULL,*ly = NULL;
1098:   PetscBool              use_gp_coords = PETSC_FALSE,set,output_gnuplot = PETSC_FALSE,glvis = PETSC_FALSE,change = PETSC_FALSE;
1099:   char                   filename[PETSC_MAX_PATH_LEN];
1100:   PetscErrorCode         ierr;


1104:   PetscOptionsGetBool(NULL,NULL,"-gnuplot",&output_gnuplot,NULL);
1105:   PetscOptionsGetBool(NULL,NULL,"-glvis",&glvis,NULL);
1106:   PetscOptionsGetBool(NULL,NULL,"-change",&change,NULL);

1108:   /* Generate the da for velocity and pressure */
1109:   /*
1110:   We use Q1 elements for the temperature.
1111:   FEM has a 9-point stencil (BOX) or connectivity pattern
1112:   Num nodes in each direction is mx+1, my+1
1113:   */
1114:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
1115:   p_dof         = P_DOFS; /* p - pressure */
1116:   dof           = u_dof+p_dof;
1117:   stencil_width = 1;
1118:   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);

1120:   DMSetMatType(da_Stokes,MATAIJ);
1121:   DMSetFromOptions(da_Stokes);
1122:   DMSetUp(da_Stokes);
1123:   DMDASetFieldName(da_Stokes,0,"Vx");
1124:   DMDASetFieldName(da_Stokes,1,"Vy");
1125:   DMDASetFieldName(da_Stokes,2,"P");

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

1130:   /* Generate element properties, we will assume all material properties are constant over the element */
1131:   /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!!  */
1132:   DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1133:   DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);

1135:   prop_dof           = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1136:   prop_stencil_width = 0;
1137:   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);
1138:   DMSetFromOptions(da_prop);
1139:   DMSetUp(da_prop);
1140:   PetscFree(lx);
1141:   PetscFree(ly);

1143:   /* define centroid positions */
1144:   DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1145:   dx   = 1.0/((PetscReal)(M));
1146:   dy   = 1.0/((PetscReal)(N));

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

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

1153:   DMCreateGlobalVector(da_prop,&properties);
1154:   DMCreateLocalVector(da_prop,&l_properties);
1155:   DMDAVecGetArray(da_prop,l_properties,&element_props);

1157:   DMGetCoordinateDM(da_prop,&prop_cda);
1158:   DMGetCoordinatesLocal(da_prop,&prop_coords);
1159:   DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);

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

1163:   DMGetCoordinateDM(da_Stokes,&vel_cda);
1164:   DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1165:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);

1167:   /* interpolate the coordinates */
1168:   for (j = sj; j < sj+ny; j++) {
1169:     for (i = si; i < si+nx; i++) {
1170:       PetscInt    ngp;
1171:       PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1172:       PetscScalar el_coords[8];

1174:       GetElementCoords(_vel_coords,i,j,el_coords);
1175:       ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

1177:       for (p = 0; p < GAUSS_POINTS; p++) {
1178:         PetscScalar gp_x,gp_y;
1179:         PetscInt    n;
1180:         PetscScalar xi_p[2],Ni_p[4];

1182:         xi_p[0] = gp_xi[p][0];
1183:         xi_p[1] = gp_xi[p][1];
1184:         ConstructQ12D_Ni(xi_p,Ni_p);

1186:         gp_x = 0.0;
1187:         gp_y = 0.0;
1188:         for (n = 0; n < NODES_PER_EL; n++) {
1189:           gp_x = gp_x+Ni_p[n]*el_coords[2*n];
1190:           gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1191:         }
1192:         element_props[j][i].gp_coords[2*p]   = gp_x;
1193:         element_props[j][i].gp_coords[2*p+1] = gp_y;
1194:       }
1195:     }
1196:   }

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

1201:   for (j = sj; j < sj+ny; j++) {
1202:     for (i = si; i < si+nx; i++) {
1203:       PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1204:       PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1205:       PetscReal coord_x,coord_y;

1207:       if (coefficient_structure == 0) {
1208:         PetscReal opts_eta0,opts_eta1,opts_xc;
1209:         PetscInt  opts_nz;

1211:         opts_eta0 = 1.0;
1212:         opts_eta1 = 1.0;
1213:         opts_xc   = 0.5;
1214:         opts_nz   = 1;

1216:         PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1217:         PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1218:         PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1219:         PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);

1221:         for (p = 0; p < GAUSS_POINTS; p++) {
1222:           coord_x = centroid_x;
1223:           coord_y = centroid_y;
1224:           if (use_gp_coords) {
1225:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1226:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1227:           }

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

1232:           element_props[j][i].fx[p] = 0.0;
1233:           element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1234:         }
1235:       } else if (coefficient_structure == 1) { /* square sinker */
1236:         PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;

1238:         opts_eta0 = 1.0;
1239:         opts_eta1 = 1.0;
1240:         opts_dx   = 0.50;
1241:         opts_dy   = 0.50;

1243:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1244:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1245:         PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1246:         PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);

1248:         for (p = 0; p < GAUSS_POINTS; p++) {
1249:           coord_x = centroid_x;
1250:           coord_y = centroid_y;
1251:           if (use_gp_coords) {
1252:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1253:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1254:           }

1256:           element_props[j][i].eta[p] = opts_eta0;
1257:           element_props[j][i].fx[p]  = 0.0;
1258:           element_props[j][i].fy[p]  = 0.0;

1260:           if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1261:             if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1262:               element_props[j][i].eta[p] =  opts_eta1;
1263:               element_props[j][i].fx[p]  =  0.0;
1264:               element_props[j][i].fy[p]  = -1.0;
1265:             }
1266:           }
1267:         }
1268:       } else if (coefficient_structure == 2) { /* circular sinker */
1269:         PetscReal opts_eta0,opts_eta1,opts_r,radius2;

1271:         opts_eta0 = 1.0;
1272:         opts_eta1 = 1.0;
1273:         opts_r    = 0.25;

1275:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1276:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1277:         PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);

1279:         for (p = 0; p < GAUSS_POINTS; p++) {
1280:           coord_x = centroid_x;
1281:           coord_y = centroid_y;
1282:           if (use_gp_coords) {
1283:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1284:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1285:           }

1287:           element_props[j][i].eta[p] = opts_eta0;
1288:           element_props[j][i].fx[p]  = 0.0;
1289:           element_props[j][i].fy[p]  = 0.0;

1291:           radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1292:           if (radius2 < opts_r*opts_r) {
1293:             element_props[j][i].eta[p] =  opts_eta1;
1294:             element_props[j][i].fx[p]  =  0.0;
1295:             element_props[j][i].fy[p]  = -1.0;
1296:           }
1297:         }
1298:       } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1299:         PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;

1301:         opts_eta0 = 1.0;
1302:         opts_eta1 = 1.0;
1303:         opts_r    = 0.25;
1304:         opts_c0x  = 0.35;       /* circle center */
1305:         opts_c0y  = 0.35;
1306:         opts_s0x  = 0.7;       /* square center */
1307:         opts_s0y  = 0.7;
1308:         opts_dx   = 0.25;
1309:         opts_dy   = 0.25;
1310:         opts_phi  = 25;

1312:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1313:         PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1314:         PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);
1315:         PetscOptionsGetReal(NULL,NULL,"-sinker_c0x",&opts_c0x,NULL);
1316:         PetscOptionsGetReal(NULL,NULL,"-sinker_c0y",&opts_c0y,NULL);
1317:         PetscOptionsGetReal(NULL,NULL,"-sinker_s0x",&opts_s0x,NULL);
1318:         PetscOptionsGetReal(NULL,NULL,"-sinker_s0y",&opts_s0y,NULL);
1319:         PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1320:         PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
1321:         PetscOptionsGetReal(NULL,NULL,"-sinker_phi",&opts_phi,NULL);
1322:         opts_phi *= PETSC_PI / 180;

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

1332:           element_props[j][i].eta[p] = opts_eta0;
1333:           element_props[j][i].fx[p]  = 0.0;
1334:           element_props[j][i].fy[p]  = -0.2;

1336:           radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1337:           if (radius2 < opts_r*opts_r
1338:               || (PetscAbs(+(coord_x - opts_s0x)*PetscCosReal(opts_phi) + (coord_y - opts_s0y)*PetscSinReal(opts_phi)) < opts_dx/2
1339:                   && PetscAbs(-(coord_x - opts_s0x)*PetscSinReal(opts_phi) + (coord_y - opts_s0y)*PetscCosReal(opts_phi)) < opts_dy/2)) {
1340:             element_props[j][i].eta[p] =  opts_eta1;
1341:             element_props[j][i].fx[p]  =  0.0;
1342:             element_props[j][i].fy[p]  = -1.0;
1343:           }
1344:         }
1345:       } else if (coefficient_structure == 4) { /* subdomain jump */
1346:         PetscReal opts_mag,opts_eta0;
1347:         PetscInt  opts_nz,px,py;
1348:         PetscBool jump;

1350:         opts_mag  = 1.0;
1351:         opts_eta0 = 1.0;
1352:         opts_nz   = 1;

1354:         PetscOptionsGetReal(NULL,NULL,"-jump_eta0",&opts_eta0,NULL);
1355:         PetscOptionsGetReal(NULL,NULL,"-jump_magnitude",&opts_mag,NULL);
1356:         PetscOptionsGetInt(NULL,NULL,"-jump_nz",&opts_nz,NULL);
1357:         DMDAGetInfo(da_Stokes,NULL,NULL,NULL,NULL,&px,&py,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1358:         if (px%2) {
1359:           jump = (PetscBool)(PetscGlobalRank%2);
1360:         } else {
1361:           jump = (PetscBool)((PetscGlobalRank/px)%2 ? PetscGlobalRank%2 : !(PetscGlobalRank%2));
1362:         }
1363:         for (p = 0; p < GAUSS_POINTS; p++) {
1364:           coord_x = centroid_x;
1365:           coord_y = centroid_y;
1366:           if (use_gp_coords) {
1367:             coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1368:             coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1369:           }

1371:           element_props[j][i].eta[p] = jump ? PetscPowReal(10.0,opts_mag) : opts_eta0;
1372:           element_props[j][i].fx[p]  = 0.0;
1373:           element_props[j][i].fy[p]  = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1374:         }
1375:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1376:     }
1377:   }
1378:   DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);

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

1382:   DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1383:   DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1384:   DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);

1386:   if (output_gnuplot) {
1387:     DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1388:     DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coefficients for Stokes eqn.","properties");
1389:   }

1391:   if (glvis) {
1392:     Vec         glv_prop,etaf;
1393:     PetscViewer view;
1394:     PetscInt    dim;
1395:     const char  *fec = {"FiniteElementCollection: L2_2D_P1"};

1397:     DMGetDimension(da_Stokes,&dim);
1398:     VecCreateSeq(PETSC_COMM_SELF,GAUSS_POINTS*mx*mx,&etaf);
1399:     PetscObjectSetName((PetscObject)etaf,"viscosity");
1400:     PetscViewerGLVisOpen(PETSC_COMM_WORLD,PETSC_VIEWER_GLVIS_SOCKET,NULL,PETSC_DECIDE,&view);
1401:     PetscViewerGLVisSetFields(view,1,&fec,&dim,glvis_extract_eta,(PetscObject*)&etaf,da_prop,NULL);
1402:     DMGetLocalVector(da_prop,&glv_prop);
1403:     DMGlobalToLocalBegin(da_prop,properties,INSERT_VALUES,glv_prop);
1404:     DMGlobalToLocalEnd(da_prop,properties,INSERT_VALUES,glv_prop);
1405:     VecSetDM(etaf,da_Stokes);
1406:     VecView(glv_prop,view);
1407:     PetscViewerDestroy(&view);
1408:     DMRestoreLocalVector(da_prop,&glv_prop);
1409:     VecDestroy(&etaf);
1410:   }

1412:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1413:   DMCreateMatrix(da_Stokes,&A);
1414:   DMCreateMatrix(da_Stokes,&B);
1415:   DMCreateGlobalVector(da_Stokes,&f);
1416:   DMCreateGlobalVector(da_Stokes,&X);

1418:   /* assemble A11 */
1419:   MatZeroEntries(A);
1420:   MatZeroEntries(B);
1421:   VecZeroEntries(f);

1423:   AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1424:   AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1425:   /* build force vector */
1426:   AssembleF_Stokes(f,da_Stokes,da_prop,properties);

1428:   DMDABCApplyFreeSlip(da_Stokes,A,f);
1429:   DMDABCApplyFreeSlip(da_Stokes,B,NULL);

1431:   /* SOLVE */
1432:   KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1433:   KSPSetOptionsPrefix(ksp_S,"stokes_");
1434:   KSPSetDM(ksp_S,da_Stokes);
1435:   KSPSetDMActive(ksp_S,PETSC_FALSE);
1436:   KSPSetOperators(ksp_S,A,B);
1437:   KSPSetFromOptions(ksp_S);
1438:   {
1439:     PC             pc;
1440:     const PetscInt ufields[] = {0,1},pfields[1] = {2};

1442:     KSPGetPC(ksp_S,&pc);
1443:     PCFieldSplitSetBlockSize(pc,3);
1444:     PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1445:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1446:   }

1448:   {
1449:     PC        pc;
1450:     PetscBool same = PETSC_FALSE;
1451:     KSPGetPC(ksp_S,&pc);
1452:     PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&same);
1453:     if (same) {
1454:       PetscBool usedivmat = PETSC_FALSE;
1455:       KSPSetOperators(ksp_S,A,A);

1457:       PetscOptionsGetBool(NULL,NULL,"-stokes_pc_bddc_use_divergence",&usedivmat,NULL);
1458:       if (usedivmat) {
1459:         IS       *fields,vel;
1460:         PetscInt i,nf;

1462:         DMCreateFieldDecomposition(da_Stokes,&nf,NULL,&fields,NULL);
1463:         ISConcatenate(PETSC_COMM_WORLD,2,fields,&vel);
1464:         MatZeroRowsIS(B,fields[2],1.0,NULL,NULL); /* we put 1.0 on the diagonal to pick the pressure average too */
1465:         MatTranspose(B,MAT_INPLACE_MATRIX,&B);
1466:         MatZeroRowsIS(B,vel,0.0,NULL,NULL);
1467:         ISDestroy(&vel);
1468:         PCBDDCSetDivergenceMat(pc,B,PETSC_FALSE,NULL);
1469:         for (i=0;i<nf;i++) {
1470:           ISDestroy(&fields[i]);
1471:         }
1472:         PetscFree(fields);
1473:       }
1474:     }
1475:   }

1477:   {
1478:     PC        pc_S;
1479:     KSP       *sub_ksp,ksp_U;
1480:     PetscInt  nsplits;
1481:     DM        da_U;
1482:     PetscBool is_pcfs;

1484:     KSPSetUp(ksp_S);
1485:     KSPGetPC(ksp_S,&pc_S);

1487:     is_pcfs = PETSC_FALSE;
1488:     PetscObjectTypeCompare((PetscObject)pc_S,PCFIELDSPLIT,&is_pcfs);

1490:     if (is_pcfs) {
1491:       PCFieldSplitGetSubKSP(pc_S,&nsplits,&sub_ksp);
1492:       ksp_U = sub_ksp[0];
1493:       PetscFree(sub_ksp);

1495:       if (nsplits == 2) {
1496:         DMDACreateCompatibleDMDA(da_Stokes,2,&da_U);

1498:         KSPSetDM(ksp_U,da_U);
1499:         KSPSetDMActive(ksp_U,PETSC_FALSE);
1500:         DMDestroy(&da_U);
1501:         if (change) {
1502:           const char opt[] = "-stokes_fieldsplit_u_pc_factor_mat_solver_type mumps";
1503:           PetscOptionsInsertString(NULL,opt);
1504:           KSPSetFromOptions(ksp_U);
1505:         }
1506:         KSPSetUp(ksp_U);
1507:       }
1508:     }
1509:   }

1511:   KSPSolve(ksp_S,f,X);

1513:   PetscOptionsGetString(NULL,NULL,"-o",filename,sizeof(filename),&set);
1514:   if (set) {
1515:     char        *ext;
1516:     PetscViewer viewer;
1517:     PetscBool   flg;
1518:     PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1519:     PetscStrrchr(filename,'.',&ext);
1520:     PetscStrcmp("vts",ext,&flg);
1521:     if (flg) {
1522:       PetscViewerSetType(viewer,PETSCVIEWERVTK);
1523:     } else {
1524:       PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1525:     }
1526:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1527:     PetscViewerFileSetName(viewer,filename);
1528:     VecView(X,viewer);
1529:     PetscViewerDestroy(&viewer);
1530:   }
1531:   if (output_gnuplot) {
1532:     DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
1533:   }

1535:   if (glvis) {
1536:     PetscViewer view;

1538:     PetscViewerCreate(PETSC_COMM_WORLD,&view);
1539:     PetscViewerSetType(view,PETSCVIEWERGLVIS);
1540:     VecView(X,view);
1541:     PetscViewerDestroy(&view);
1542:   }

1544:   KSPGetIterationNumber(ksp_S,&its);

1546:   if (coefficient_structure == 0) {
1547:     PetscReal opts_eta0,opts_eta1,opts_xc;
1548:     PetscInt  opts_nz,N;
1549:     DM        da_Stokes_analytic;
1550:     Vec       X_analytic;
1551:     PetscReal nrm1[3],nrm2[3],nrmI[3];

1553:     opts_eta0 = 1.0;
1554:     opts_eta1 = 1.0;
1555:     opts_xc   = 0.5;
1556:     opts_nz   = 1;

1558:     PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1559:     PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1560:     PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1561:     PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);

1563:     DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1564:     if (output_gnuplot) {
1565:       DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
1566:     }
1567:     DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);

1569:     VecAXPY(X_analytic,-1.0,X);
1570:     VecGetSize(X_analytic,&N);
1571:     N    = N/3;

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

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

1581:     VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1582:     VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1583:     VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);

1585:     DMDestroy(&da_Stokes_analytic);
1586:     VecDestroy(&X_analytic);
1587:   }

1589:   KSPDestroy(&ksp_S);
1590:   VecDestroy(&X);
1591:   VecDestroy(&f);
1592:   MatDestroy(&A);
1593:   MatDestroy(&B);

1595:   DMDestroy(&da_Stokes);
1596:   DMDestroy(&da_prop);

1598:   VecDestroy(&properties);
1599:   VecDestroy(&l_properties);
1600:   return(0);
1601: }

1603: int main(int argc,char **args)
1604: {
1606:   PetscInt       mx,my;

1608:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
1609:   mx   = my = 10;
1610:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1611:   PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1612:   solve_stokes_2d_coupled(mx,my);
1613:   PetscFinalize();
1614:   return ierr;
1615: }

1617: /* -------------------------- helpers for boundary conditions -------------------------------- */
1618: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
1619: {
1620:   DM                     cda;
1621:   Vec                    coords;
1622:   PetscInt               si,sj,nx,ny,i,j;
1623:   PetscInt               M,N;
1624:   DMDACoor2d             **_coords;
1625:   const PetscInt         *g_idx;
1626:   PetscInt               *bc_global_ids;
1627:   PetscScalar            *bc_vals;
1628:   PetscInt               nbcs;
1629:   PetscInt               n_dofs;
1630:   PetscErrorCode         ierr;
1631:   ISLocalToGlobalMapping ltogm;

1634:   DMGetLocalToGlobalMapping(da,&ltogm);
1635:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1637:   DMGetCoordinateDM(da,&cda);
1638:   DMGetCoordinatesLocal(da,&coords);
1639:   DMDAVecGetArray(cda,coords,&_coords);
1640:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1641:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1643:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
1644:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

1649:   i = nx-1;
1650:   for (j = 0; j < ny; j++) {
1651:     PetscInt local_id;

1653:     local_id = i+j*nx;

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

1657:     bc_vals[j] =  0.0;
1658:   }
1659:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1660:   nbcs = 0;
1661:   if ((si+nx) == (M)) nbcs = ny;

1663:   if (b) {
1664:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1665:     VecAssemblyBegin(b);
1666:     VecAssemblyEnd(b);
1667:   }
1668:   if (A) {
1669:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1670:   }

1672:   PetscFree(bc_vals);
1673:   PetscFree(bc_global_ids);

1675:   DMDAVecRestoreArray(cda,coords,&_coords);
1676:   return(0);
1677: }

1679: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
1680: {
1681:   DM                     cda;
1682:   Vec                    coords;
1683:   PetscInt               si,sj,nx,ny,i,j;
1684:   PetscInt               M,N;
1685:   DMDACoor2d             **_coords;
1686:   const PetscInt         *g_idx;
1687:   PetscInt               *bc_global_ids;
1688:   PetscScalar            *bc_vals;
1689:   PetscInt               nbcs;
1690:   PetscInt               n_dofs;
1691:   PetscErrorCode         ierr;
1692:   ISLocalToGlobalMapping ltogm;

1695:   DMGetLocalToGlobalMapping(da,&ltogm);
1696:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1698:   DMGetCoordinateDM(da,&cda);
1699:   DMGetCoordinatesLocal(da,&coords);
1700:   DMDAVecGetArray(cda,coords,&_coords);
1701:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1702:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1704:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
1705:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

1710:   i = 0;
1711:   for (j = 0; j < ny; j++) {
1712:     PetscInt local_id;

1714:     local_id = i+j*nx;

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

1718:     bc_vals[j] =  0.0;
1719:   }
1720:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1721:   nbcs = 0;
1722:   if (si == 0) nbcs = ny;

1724:   if (b) {
1725:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1726:     VecAssemblyBegin(b);
1727:     VecAssemblyEnd(b);
1728:   }

1730:   if (A) {
1731:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1732:   }

1734:   PetscFree(bc_vals);
1735:   PetscFree(bc_global_ids);

1737:   DMDAVecRestoreArray(cda,coords,&_coords);
1738:   return(0);
1739: }

1741: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
1742: {
1743:   DM                     cda;
1744:   Vec                    coords;
1745:   PetscInt               si,sj,nx,ny,i,j;
1746:   PetscInt               M,N;
1747:   DMDACoor2d             **_coords;
1748:   const PetscInt         *g_idx;
1749:   PetscInt               *bc_global_ids;
1750:   PetscScalar            *bc_vals;
1751:   PetscInt               nbcs;
1752:   PetscInt               n_dofs;
1753:   PetscErrorCode         ierr;
1754:   ISLocalToGlobalMapping ltogm;

1757:   DMGetLocalToGlobalMapping(da,&ltogm);
1758:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1760:   DMGetCoordinateDM(da,&cda);
1761:   DMGetCoordinatesLocal(da,&coords);
1762:   DMDAVecGetArray(cda,coords,&_coords);
1763:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1764:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1766:   PetscMalloc1(nx,&bc_global_ids);
1767:   PetscMalloc1(nx,&bc_vals);

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

1772:   j = ny-1;
1773:   for (i = 0; i < nx; i++) {
1774:     PetscInt local_id;

1776:     local_id = i+j*nx;

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

1780:     bc_vals[i] =  0.0;
1781:   }
1782:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1783:   nbcs = 0;
1784:   if ((sj+ny) == (N)) nbcs = nx;

1786:   if (b) {
1787:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1788:     VecAssemblyBegin(b);
1789:     VecAssemblyEnd(b);
1790:   }
1791:   if (A) {
1792:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1793:   }

1795:   PetscFree(bc_vals);
1796:   PetscFree(bc_global_ids);

1798:   DMDAVecRestoreArray(cda,coords,&_coords);
1799:   return(0);
1800: }

1802: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
1803: {
1804:   DM                     cda;
1805:   Vec                    coords;
1806:   PetscInt               si,sj,nx,ny,i,j;
1807:   PetscInt               M,N;
1808:   DMDACoor2d             **_coords;
1809:   const PetscInt         *g_idx;
1810:   PetscInt               *bc_global_ids;
1811:   PetscScalar            *bc_vals;
1812:   PetscInt               nbcs;
1813:   PetscInt               n_dofs;
1814:   PetscErrorCode         ierr;
1815:   ISLocalToGlobalMapping ltogm;

1818:   DMGetLocalToGlobalMapping(da,&ltogm);
1819:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1821:   DMGetCoordinateDM(da,&cda);
1822:   DMGetCoordinatesLocal(da,&coords);
1823:   DMDAVecGetArray(cda,coords,&_coords);
1824:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1825:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1827:   PetscMalloc1(nx,&bc_global_ids);
1828:   PetscMalloc1(nx,&bc_vals);

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

1833:   j = 0;
1834:   for (i = 0; i < nx; i++) {
1835:     PetscInt local_id;

1837:     local_id = i+j*nx;

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

1841:     bc_vals[i] =  0.0;
1842:   }
1843:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1844:   nbcs = 0;
1845:   if (sj == 0) nbcs = nx;

1847:   if (b) {
1848:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1849:     VecAssemblyBegin(b);
1850:     VecAssemblyEnd(b);
1851:   }
1852:   if (A) {
1853:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1854:   }

1856:   PetscFree(bc_vals);
1857:   PetscFree(bc_global_ids);

1859:   DMDAVecRestoreArray(cda,coords,&_coords);
1860:   return(0);
1861: }

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

1869:   BCApplyZero_NORTH(da_Stokes,1,A,f);
1870:   BCApplyZero_EAST(da_Stokes,0,A,f);
1871:   BCApplyZero_SOUTH(da_Stokes,1,A,f);
1872:   BCApplyZero_WEST(da_Stokes,0,A,f);
1873:   return(0);
1874: }

1876: /*TEST

1878:    build:
1879:       requires: !complex !single

1881:    test:
1882:       args: -stokes_pc_use_amat -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type lu -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short

1884:    testset:
1885:       args: -stokes_pc_use_amat -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_fieldsplit_u_ksp_type preonly -stokes_fieldsplit_u_pc_type lu -stokes_fieldsplit_p_ksp_type preonly -stokes_fieldsplit_p_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short
1886:       test:
1887:          suffix: 2
1888:          args:
1889:          output_file: output/ex43_1.out
1890:       test:
1891:          requires: mumps
1892:          suffix: 2_mumps
1893:          args: -change true -stokes_ksp_view
1894:          output_file: output/ex43_2_mumps.out
1895:          # mumps INFO,INFOG,RINFO,RINFOG may vary on different archs, so keep just a stable one
1896:          filter:  egrep -v "(INFOG\([^7]|INFO\(|\[0\])"

1898:    test:
1899:       suffix: 3
1900:       nsize: 4
1901:       args: -stokes_pc_use_amat -stokes_ksp_type gcr -stokes_ksp_gcr_restart 60 -stokes_ksp_norm_type unpreconditioned -stokes_ksp_rtol 1.e-2 -c_str 3 -sinker_eta0 1.0 -sinker_eta1 100 -sinker_dx 0.4 -sinker_dy 0.3 -mx 128 -my 128 -stokes_ksp_monitor_short -stokes_pc_type mg -stokes_mg_levels_pc_type fieldsplit -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_levels_pc_fieldsplit_block_size 3 -stokes_mg_levels_pc_fieldsplit_0_fields 0,1 -stokes_mg_levels_pc_fieldsplit_1_fields 2 -stokes_mg_levels_fieldsplit_0_pc_type sor -stokes_mg_levels_fieldsplit_1_pc_type sor -stokes_mg_levels_ksp_type chebyshev -stokes_mg_levels_ksp_max_it 1 -stokes_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -stokes_pc_mg_levels 4 -stokes_ksp_view

1903:    test:
1904:       suffix: 4
1905:       nsize: 4
1906:       args: -stokes_ksp_type pipegcr -stokes_ksp_pipegcr_mmax 60 -stokes_ksp_pipegcr_unroll_w 1 -stokes_ksp_norm_type natural -c_str 3 -sinker_eta0 1.0 -sinker_eta1 100 -sinker_dx 0.4 -sinker_dy 0.3 -mx 128 -my 128 -stokes_ksp_monitor_short -stokes_pc_type mg -stokes_mg_levels_pc_type fieldsplit -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_levels_pc_fieldsplit_block_size 3 -stokes_mg_levels_pc_fieldsplit_0_fields 0,1 -stokes_mg_levels_pc_fieldsplit_1_fields 2 -stokes_mg_levels_fieldsplit_0_pc_type sor -stokes_mg_levels_fieldsplit_1_pc_type sor -stokes_mg_levels_ksp_type chebyshev -stokes_mg_levels_ksp_max_it 1 -stokes_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -stokes_pc_mg_levels 4 -stokes_ksp_view

1908:    test:
1909:       suffix: 5
1910:       nsize: 4
1911:       args: -stokes_pc_fieldsplit_off_diag_use_amat -stokes_ksp_type pipegcr -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type bjacobi -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type bjacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short -stokes_ksp_view

1913:    test:
1914:       suffix: 6
1915:       nsize: 8
1916:       args: -stokes_ksp_view -stokes_pc_type mg -stokes_pc_mg_levels 2 -stokes_mg_coarse_pc_type telescope -stokes_mg_coarse_pc_telescope_reduction_factor 2 -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_coarse_pc_telescope_subcomm_type contiguous

1918:    test:
1919:       suffix: bjacobi
1920:       nsize: 4
1921:       args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type aij -stokes_ksp_converged_reason

1923:    test:
1924:       suffix: bjacobi_baij
1925:       nsize: 4
1926:       args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type baij -stokes_ksp_converged_reason
1927:       output_file: output/ex43_bjacobi.out

1929:    test:
1930:       suffix: nested_gmg
1931:       nsize: 4
1932:       args: -stokes_pc_fieldsplit_off_diag_use_amat -mx 16 -my 16 -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_fieldsplit_u_pc_type mg -stokes_fieldsplit_u_pc_mg_levels 5 -stokes_fieldsplit_u_pc_mg_galerkin pmat -stokes_fieldsplit_u_ksp_type cg -stokes_fieldsplit_u_ksp_rtol 1.0e-4 -stokes_fieldsplit_u_mg_levels_pc_type jacobi -solcx_eta0 1.0e4 -stokes_fieldsplit_u_ksp_converged_reason -stokes_ksp_converged_reason -stokes_fieldsplit_p_sub_pc_factor_zeropivot 1.e-8

1934:    test:
1935:       suffix: fetidp
1936:       nsize: 8
1937:       args: -dm_mat_type is -stokes_ksp_type fetidp -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_converged_reason -stokes_fetidp_pc_fieldsplit_schur_fact_type diag -stokes_fetidp_fieldsplit_p_pc_type bjacobi -stokes_fetidp_fieldsplit_lag_ksp_type preonly -stokes_fetidp_fieldsplit_p_ksp_type preonly -stokes_ksp_fetidp_pressure_field 2 -stokes_fetidp_pc_fieldsplit_schur_scale -1

1939:    test:
1940:       suffix: fetidp_unsym
1941:       nsize: 8
1942:       args: -dm_mat_type is -stokes_ksp_type fetidp -stokes_ksp_monitor_true_residual -stokes_ksp_converged_reason -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd

1944:    test:
1945:       suffix: bddc_stokes_deluxe
1946:       nsize: 8
1947:       args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_coarse_redundant_pc_type svd -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc

1949:    test:
1950:       suffix: bddc_stokes_subdomainjump_deluxe
1951:       nsize: 9
1952:       args: -c_str 4 -jump_magnitude 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_coarse_redundant_pc_type svd -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc

1954: TEST*/