Actual source code: ex43.c

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

 30: /* Contributed by Dave May */

 32: #include <petscksp.h>
 33: #include <petscdmda.h>

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

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


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

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

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

 60: typedef struct {
 61:   PetscScalar u_dof;
 62:   PetscScalar v_dof;
 63:   PetscScalar p_dof;
 64: } StokesDOF;


 67: /*

 69: D = [ 2.eta   0   0   ]
 70: [   0   2.eta 0   ]
 71: [   0     0   eta ]

 73: B = [ d_dx   0   ]
 74: [  0    d_dy ]
 75: [ d_dy  d_dx ]

 77: */

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

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

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

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

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

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

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

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

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


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

143:   if (det_J != NULL) {
144:     *det_J = J;
145:   }
146: }

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


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

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

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

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

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

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

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

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

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

225:   return(0);
226: }

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

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

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

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


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

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

280:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

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

284:   proc_J = rank/cpu_x;
285:   proc_I = rank-cpu_x*proc_J;

287:   PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);
288:   PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);

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

295:   VecCreate(PETSC_COMM_WORLD,&vly);
296:   VecSetSizes(vly,PETSC_DECIDE,cpu_y);
297:   VecSetFromOptions(vly);

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



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

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



330:   *_lx = LX;
331:   *_ly = LY;

333:   VecDestroy(&vlx);
334:   VecDestroy(&vly);

336:   return(0);
337: }

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

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,"### Element geometry for processor %1.4d ### \n",rank);

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

375:   PetscFClose(PETSC_COMM_SELF,fp);
376:   return(0);
377: }

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

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

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


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

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


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

427:       coord_x = _coords[j][i].x;
428:       coord_y = _coords[j][i].y;

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

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

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

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

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

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


478:   DMDAGetCoordinateDA(da,&cda);
479:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);

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


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

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

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

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

505:   PetscFClose(PETSC_COMM_SELF,fp);
506:   return(0);
507: }


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

517:   nc = u_NPE*u_dof;

519:   r = w_dof*wi+wd;
520:   c = u_dof*ui+ud;

522:   ij = r*nc+c;

524:   return ij;
525: }

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


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

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

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

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


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

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

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


587:   /* define quadrature rule */
588:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

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

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

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

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

617:   PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
618:   FormGradientOperatorQ1(Ge,coords);

620:   nr_g = U_DOFS*NODES_PER_EL;
621:   nc_g = P_DOFS*NODES_PER_EL;

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

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


641:   /* define quadrature rule */
642:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

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

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

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

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


683:   /* define quadrature rule */
684:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

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

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

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

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


725:   /* define quadrature rule */
726:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

908:   for (n = 0; n < 4; n++) {
909:     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  ];
910:     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];
911:     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    ];
912:   }
913:   return(0);
914: }

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

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

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

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


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

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

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


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

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

976:       DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
977:     }
978:   }

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


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

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

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

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


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


1018:   DMDAGetGhostedCoordinates(da,&coords);
1019:   DMDAGetCoordinateDA(da,&cda);
1020:   DMDAVecGetArray(cda,coords,&_coords);

1022:   DMCreateGlobalVector(da,&X);
1023:   DMCreateLocalVector(da,&local_X);
1024:   DMDAVecGetArray(da,local_X,&_stokes);

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

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

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

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

1044:   DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1045:   DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);

1047:   VecDestroy(&local_X);

1049:   *_da = da;
1050:   *_X  = X;

1052:   return(0);
1053: }

1057: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
1058: {
1060:   /* get the nodal fields */
1061:   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;
1062:   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;
1063:   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;
1064:   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;
1065:   return(0);
1066: }

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

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

1094:   /* define quadrature rule */
1095:   ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

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

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

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

1114:   DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1115:   DMDAGetBoundingBox(stokes_da,xymin,xymax);

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

1119:   tp_L2 = tu_L2 = tu_H1 = 0.0;

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

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

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

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

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

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

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

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


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

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

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

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

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


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

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

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

1240:   /* define centroid positions */
1241:   DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1242:   dx   = 1.0/((PetscReal)(M));
1243:   dy   = 1.0/((PetscReal)(N));

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

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

1251:   DMCreateGlobalVector(da_prop,&properties);
1252:   DMCreateLocalVector(da_prop,&l_properties);
1253:   DMDAVecGetArray(da_prop,l_properties,&element_props);

1255:   DMDAGetCoordinateDA(da_prop,&prop_cda);
1256:   DMDAGetGhostedCoordinates(da_prop,&prop_coords);
1257:   DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);

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

1261:   DMDAGetCoordinateDA(da_Stokes,&vel_cda);
1262:   DMDAGetGhostedCoordinates(da_Stokes,&vel_coords);
1263:   DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);


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

1273:       GetElementCoords(_vel_coords,i,j,el_coords);
1274:       ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);

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

1281:         xi_p[0] = gp_xi[p][0];
1282:         xi_p[1] = gp_xi[p][1];
1283:         ConstructQ12D_Ni(xi_p,Ni_p);

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

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

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

1306:       if (coefficient_structure == 0) {
1307:         PetscReal opts_eta0,opts_eta1,opts_xc;
1308:         PetscInt  opts_nz;

1310:         opts_eta0 = 1.0;
1311:         opts_eta1 = 1.0;
1312:         opts_xc   = 0.5;
1313:         opts_nz   = 1;

1315:         PetscOptionsGetReal(PETSC_NULL,"-solcx_eta0",&opts_eta0,0);
1316:         PetscOptionsGetReal(PETSC_NULL,"-solcx_eta1",&opts_eta1,0);
1317:         PetscOptionsGetReal(PETSC_NULL,"-solcx_xc",&opts_xc,0);
1318:         PetscOptionsGetInt(PETSC_NULL,"-solcx_nz",&opts_nz,0);

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


1329:           element_props[j][i].eta[p] = opts_eta0;
1330:           if (coord_x > opts_xc) {
1331:             element_props[j][i].eta[p] = opts_eta1;
1332:           }

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

1340:         opts_eta0 = 1.0;
1341:         opts_eta1 = 1.0;
1342:         opts_dx   = 0.50;
1343:         opts_dy   = 0.50;

1345:         PetscOptionsGetReal(PETSC_NULL,"-sinker_eta0",&opts_eta0,0);
1346:         PetscOptionsGetReal(PETSC_NULL,"-sinker_eta1",&opts_eta1,0);
1347:         PetscOptionsGetReal(PETSC_NULL,"-sinker_dx",&opts_dx,0);
1348:         PetscOptionsGetReal(PETSC_NULL,"-sinker_dy",&opts_dy,0);


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

1359:           element_props[j][i].eta[p] = opts_eta0;
1360:           element_props[j][i].fx[p]  = 0.0;
1361:           element_props[j][i].fy[p]  = 0.0;

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

1374:         opts_eta0 = 1.0;
1375:         opts_eta1 = 1.0;
1376:         opts_r    = 0.25;

1378:         PetscOptionsGetReal(PETSC_NULL,"-sinker_eta0",&opts_eta0,0);
1379:         PetscOptionsGetReal(PETSC_NULL,"-sinker_eta1",&opts_eta1,0);
1380:         PetscOptionsGetReal(PETSC_NULL,"-sinker_r",&opts_r,0);

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

1390:           element_props[j][i].eta[p] = opts_eta0;
1391:           element_props[j][i].fx[p]  = 0.0;
1392:           element_props[j][i].fy[p]  = 0.0;

1394:           radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1395:           if (radius2 < opts_r*opts_r) {
1396:             element_props[j][i].eta[p] =  opts_eta1;
1397:             element_props[j][i].fx[p]  =  0.0;
1398:             element_props[j][i].fy[p]  = -1.0;
1399:           }
1400:         }
1401:       } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1402:     }
1403:   }
1404:   DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);

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

1408:   DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1409:   DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1410:   DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);


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


1417:   /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1418:   DMGetMatrix(da_Stokes,MATAIJ,&A);
1419:   DMGetMatrix(da_Stokes,MATAIJ,&B);
1420:   MatGetVecs(A,&f,&X);

1422:   /* assemble A11 */
1423:   MatZeroEntries(A);
1424:   MatZeroEntries(B);
1425:   VecZeroEntries(f);

1427:   AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1428:   AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1429:   /* build force vector */
1430:   AssembleF_Stokes(f,da_Stokes,da_prop,properties);

1432:   DMDABCApplyFreeSlip(da_Stokes,A,f);
1433:   DMDABCApplyFreeSlip(da_Stokes,B,PETSC_NULL);

1435:   /* SOLVE */
1436:   KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1437:   KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
1438:   KSPSetOperators(ksp_S,A,B,SAME_NONZERO_PATTERN);
1439:   KSPSetFromOptions(ksp_S);
1440:   {
1441:     PC pc;
1442:     const PetscInt ufields[] = {0,1},pfields[1] = {2};
1443:     KSPGetPC(ksp_S,&pc);
1444:     PCFieldSplitSetFields(pc,"u",2,ufields);
1445:     PCFieldSplitSetFields(pc,"p",1,pfields);
1446:   }

1448:   KSPSolve(ksp_S,f,X);
1449:   DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");

1451:   KSPGetIterationNumber(ksp_S,&its);

1453:   /*
1454:   {
1455:   Vec x;
1456:   PetscInt L;
1457:   VecDuplicate(f,&x);
1458:   MatMult(A,X, x );
1459:   VecAXPY( x, -1.0, f );
1460:   VecNorm( x, NORM_2, &nrm );
1461:   PetscPrintf( PETSC_COMM_WORLD, "Its. %1.4d, norm |AX-f| = %1.5e \n", its, nrm );
1462:   VecDestroy(&x);

1464:   VecNorm( X, NORM_2, &nrm );
1465:   VecGetSize( X, &L );
1466:   PetscPrintf( PETSC_COMM_WORLD, "           norm |X|/sqrt{N} = %1.5e \n", nrm/sqrt( (PetscScalar)L ) );
1467:   }
1468:   */

1470:   if (coefficient_structure == 0) {
1471:     PetscReal   opts_eta0,opts_eta1,opts_xc;
1472:     PetscInt    opts_nz,N;
1473:     DM          da_Stokes_analytic;
1474:     Vec         X_analytic;
1475:     PetscReal   nrm1[3],nrm2[3],nrmI[3];

1477:     opts_eta0 = 1.0;
1478:     opts_eta1 = 1.0;
1479:     opts_xc   = 0.5;
1480:     opts_nz   = 1;

1482:     PetscOptionsGetReal(PETSC_NULL,"-solcx_eta0",&opts_eta0,0);
1483:     PetscOptionsGetReal(PETSC_NULL,"-solcx_eta1",&opts_eta1,0);
1484:     PetscOptionsGetReal(PETSC_NULL,"-solcx_xc",&opts_xc,0);
1485:     PetscOptionsGetInt(PETSC_NULL,"-solcx_nz",&opts_nz,0);


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

1491:     DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);


1494:     VecAXPY(X_analytic,-1.0,X);
1495:     VecGetSize(X_analytic,&N);
1496:     N    = N/3;

1498:     VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1499:     VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1500:     VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);

1502:     VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1503:     VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1504:     VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);

1506:     VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1507:     VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1508:     VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);

1510:     /*
1511:     PetscPrintf( PETSC_COMM_WORLD, "%1.4e   %1.4e %1.4e %1.4e    %1.4e %1.4e %1.4e    %1.4e %1.4e %1.4e\n",
1512:     1.0/mx,
1513:     nrm1[0]/(double)N, sqrt(nrm2[0]/(double)N), nrmI[0],
1514:     nrm1[1]/(double)N, sqrt(nrm2[1]/(double)N), nrmI[1],
1515:     nrm1[2]/(double)N, sqrt(nrm2[2]/(double)N), nrmI[2] );
1516:     */
1517:     DMDestroy(&da_Stokes_analytic);
1518:     VecDestroy(&X_analytic);
1519:   }


1522:   KSPDestroy(&ksp_S);
1523:   VecDestroy(&X);
1524:   VecDestroy(&f);
1525:   MatDestroy(&A);
1526:   MatDestroy(&B);

1528:   DMDestroy(&da_Stokes);
1529:   DMDestroy(&da_prop);

1531:   VecDestroy(&properties);
1532:   VecDestroy(&l_properties);

1534:   return(0);
1535: }

1539: int main(int argc,char **args)
1540: {
1542:   PetscInt       mx,my;

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

1546:   mx   = my = 10;
1547:   PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,PETSC_NULL);
1548:   PetscOptionsGetInt(PETSC_NULL,"-my",&my,PETSC_NULL);

1550:   solve_stokes_2d_coupled(mx,my);

1552:   PetscFinalize();
1553:   return 0;
1554: }

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

1560: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1561: {
1562:   DM             cda;
1563:   Vec            coords;
1564:   PetscInt       si,sj,nx,ny,i,j;
1565:   PetscInt       M,N;
1566:   DMDACoor2d       **_coords;
1567:   PetscInt       *g_idx;
1568:   PetscInt       *bc_global_ids;
1569:   PetscScalar    *bc_vals;
1570:   PetscInt       nbcs;
1571:   PetscInt       n_dofs;

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

1578:   DMDAGetCoordinateDA(da,&cda);
1579:   DMDAGetGhostedCoordinates(da,&coords);
1580:   DMDAVecGetArray(cda,coords,&_coords);
1581:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1582:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1584:   /* /// */

1586:   PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1587:   PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);

1589:   /* init the entries to -1 so VecSetValues will ignore them */
1590:   for (i = 0; i < ny*n_dofs; i++) {
1591:     bc_global_ids[i] = -1;
1592:   }

1594:   i = nx-1;
1595:   for (j = 0; j < ny; j++) {
1596:     PetscInt    local_id;

1598:     local_id = i+j*nx;

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

1602:     bc_vals[j] =  bc_val;
1603:   }
1604:   nbcs = 0;
1605:   if ((si+nx) == (M)) {
1606:     nbcs = ny;
1607:   }

1609:   if (b != PETSC_NULL) {
1610:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1611:     VecAssemblyBegin(b);
1612:     VecAssemblyEnd(b);
1613:   }
1614:   if (A != PETSC_NULL) {
1615:     MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1616:   }


1619:   PetscFree(bc_vals);
1620:   PetscFree(bc_global_ids);

1622:   DMDAVecRestoreArray(cda,coords,&_coords);
1623:   return(0);
1624: }

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

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

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

1652:   /* /// */

1654:   PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1655:   PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);

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

1662:   i = 0;
1663:   for (j = 0; j < ny; j++) {
1664:     PetscInt    local_id;

1666:     local_id = i+j*nx;

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

1670:     bc_vals[j] =  bc_val;
1671:   }
1672:   nbcs = 0;
1673:   if (si == 0) {
1674:     nbcs = ny;
1675:   }

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


1687:   PetscFree(bc_vals);
1688:   PetscFree(bc_global_ids);

1690:   DMDAVecRestoreArray(cda,coords,&_coords);
1691:   return(0);
1692: }

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

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

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

1720:   /* /// */

1722:   PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1723:   PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);

1725:   /* init the entries to -1 so VecSetValues will ignore them */
1726:   for (i = 0; i < nx; i++) {
1727:     bc_global_ids[i] = -1;
1728:   }

1730:   j = ny-1;
1731:   for (i = 0; i < nx; i++) {
1732:     PetscInt    local_id;

1734:     local_id = i+j*nx;

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

1738:     bc_vals[i] =  bc_val;
1739:   }
1740:   nbcs = 0;
1741:   if ((sj+ny) == (N)) {
1742:     nbcs = nx;
1743:   }

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


1755:   PetscFree(bc_vals);
1756:   PetscFree(bc_global_ids);

1758:   DMDAVecRestoreArray(cda,coords,&_coords);
1759:   return(0);
1760: }

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

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

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

1788:   /* /// */

1790:   PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1791:   PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);

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

1798:   j = 0;
1799:   for (i = 0; i < nx; i++) {
1800:     PetscInt    local_id;

1802:     local_id = i+j*nx;

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

1806:     bc_vals[i] =  bc_val;
1807:   }
1808:   nbcs = 0;
1809:   if (sj == 0) {
1810:     nbcs = nx;
1811:   }


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


1824:   PetscFree(bc_vals);
1825:   PetscFree(bc_global_ids);

1827:   DMDAVecRestoreArray(cda,coords,&_coords);
1828:   return(0);
1829: }

1833: /*
1834: Free slip sides.
1835: */
1838: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1839: {

1843:   BCApply_NORTH(da_Stokes,1,0.0,A,f);
1844:   BCApply_EAST(da_Stokes,0,0.0,A,f);
1845:   BCApply_SOUTH(da_Stokes,1,0.0,A,f);
1846:   BCApply_WEST(da_Stokes,0,0.0,A,f);
1847:   return(0);
1848: }