Actual source code: ex43.c
petsc-3.5.4 2015-05-23
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: -mx : number elements in x-direciton \n\
6: -my : number elements in y-direciton \n\
7: -c_str : indicates the structure of the coefficients to use. \n\
8: -c_str 0 => Setup for an analytic solution with a vertical jump in viscosity. This problem is driven by the \n\
9: forcing function f = (0, sin(n_z pi y)cos(pi x). \n\
10: Parameters: \n\
11: -solcx_eta0 : the viscosity to the left of the interface \n\
12: -solcx_eta1 : the viscosity to the right of the interface \n\
13: -solcx_xc : the location of the interface \n\
14: -solcx_nz : the wavenumber in the y direction \n\
15: -c_str 1 => Setup for a rectangular blob located in the origin of the domain. \n\
16: Parameters: \n\
17: -sinker_eta0 : the viscosity of the background fluid \n\
18: -sinker_eta1 : the viscosity of the blob \n\
19: -sinker_dx : the width of the blob \n\
20: -sinker_dy : the width of the blob \n\
21: -c_str 2 => Setup for a circular blob located in the origin of the domain. \n\
22: Parameters: \n\
23: -sinker_eta0 : the viscosity of the background fluid \n\
24: -sinker_eta1 : the viscosity of the blob \n\
25: -sinker_r : radius of the blob \n\
26: -c_str 3 => Circular and rectangular inclusion\n\
27: Parameters as for cases 1 and 2 (above)\n\
28: -use_gp_coords : evaluate the viscosity and the body force at the global coordinates of the quadrature points.\n\
29: By default, eta and the body force are evaulated at the element center and applied as a constant over the entire element.\n";
31: /* Contributed by Dave May */
33: #include <petscksp.h>
34: #include <petscdm.h>
35: #include <petscdmda.h>
37: /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
38: #include "ex43-solcx.h"
40: static PetscErrorCode DMDABCApplyFreeSlip(DM,Mat,Vec);
43: #define NSD 2 /* number of spatial dimensions */
44: #define NODES_PER_EL 4 /* nodes per element */
45: #define U_DOFS 2 /* degrees of freedom per velocity node */
46: #define P_DOFS 1 /* degrees of freedom per pressure node */
47: #define GAUSS_POINTS 4
49: /* cell based evaluation */
50: typedef struct {
51: PetscScalar eta,fx,fy;
52: } Coefficients;
54: /* Gauss point based evaluation 8+4+4+4 = 20 */
55: typedef struct {
56: PetscScalar gp_coords[2*GAUSS_POINTS];
57: PetscScalar eta[GAUSS_POINTS];
58: PetscScalar fx[GAUSS_POINTS];
59: PetscScalar fy[GAUSS_POINTS];
60: } GaussPointCoefficients;
62: typedef struct {
63: PetscScalar u_dof;
64: PetscScalar v_dof;
65: PetscScalar p_dof;
66: } StokesDOF;
69: /*
71: D = [ 2.eta 0 0 ]
72: [ 0 2.eta 0 ]
73: [ 0 0 eta ]
75: B = [ d_dx 0 ]
76: [ 0 d_dy ]
77: [ d_dy d_dx ]
79: */
81: /* FEM routines */
82: /*
83: Element: Local basis function ordering
84: 1-----2
85: | |
86: | |
87: 0-----3
88: */
89: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
90: {
91: PetscScalar xi = _xi[0];
92: PetscScalar eta = _xi[1];
94: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
95: Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
96: Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
97: Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
98: }
100: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
101: {
102: PetscScalar xi = _xi[0];
103: PetscScalar eta = _xi[1];
105: GNi[0][0] = -0.25*(1.0-eta);
106: GNi[0][1] = -0.25*(1.0+eta);
107: GNi[0][2] = 0.25*(1.0+eta);
108: GNi[0][3] = 0.25*(1.0-eta);
110: GNi[1][0] = -0.25*(1.0-xi);
111: GNi[1][1] = 0.25*(1.0-xi);
112: GNi[1][2] = 0.25*(1.0+xi);
113: GNi[1][3] = -0.25*(1.0+xi);
114: }
116: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
117: {
118: PetscScalar J00,J01,J10,J11,J;
119: PetscScalar iJ00,iJ01,iJ10,iJ11;
120: PetscInt i;
122: J00 = J01 = J10 = J11 = 0.0;
123: for (i = 0; i < NODES_PER_EL; i++) {
124: PetscScalar cx = coords[2*i+0];
125: PetscScalar cy = coords[2*i+1];
127: J00 = J00+GNi[0][i]*cx; /* J_xx = dx/dxi */
128: J01 = J01+GNi[0][i]*cy; /* J_xy = dy/dxi */
129: J10 = J10+GNi[1][i]*cx; /* J_yx = dx/deta */
130: J11 = J11+GNi[1][i]*cy; /* J_yy = dy/deta */
131: }
132: J = (J00*J11)-(J01*J10);
134: iJ00 = J11/J;
135: iJ01 = -J01/J;
136: iJ10 = -J10/J;
137: iJ11 = J00/J;
140: for (i = 0; i < NODES_PER_EL; i++) {
141: GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
142: GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
143: }
145: if (det_J != NULL) *det_J = J;
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 != NULL) {
175: *mxl = m;
176: if ((sx+m) == M) *mxl = m-1; /* last proc */
177: }
178: if (myl != NULL) {
179: *myl = n;
180: if ((sy+n) == N) *myl = n-1; /* last proc */
181: }
182: if (mzl != NULL) {
183: *mzl = p;
184: if ((sz+p) == P) *mzl = p-1; /* last proc */
185: }
186: return(0);
187: }
191: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
192: {
193: PetscInt si,sj,sk;
196: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
198: *sx = si;
199: if (si) *sx = si+1;
201: *sy = sj;
202: if (sj) *sy = sj+1;
204: if (sk) {
205: *sz = sk;
206: if (sk != 0) *sz = sk+1;
207: }
209: DMDAGetLocalElementSize(da,mx,my,mz);
210: return(0);
211: }
213: /*
214: i,j are the element indices
215: The unknown is a vector quantity.
216: The s[].c is used to indicate the degree of freedom.
217: */
220: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
221: {
223: /* velocity */
224: /* node 0 */
225: s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0; /* Vx0 */
226: s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1; /* Vy0 */
228: /* node 1 */
229: s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0; /* Vx1 */
230: s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1; /* Vy1 */
232: /* node 2 */
233: s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0; /* Vx2 */
234: s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1; /* Vy2 */
236: /* node 3 */
237: s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0; /* Vx3 */
238: s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1; /* Vy3 */
241: /* pressure */
242: s_p[0].i = i;s_p[0].j = j;s_p[0].c = 2; /* P0 */
243: s_p[1].i = i;s_p[1].j = j+1;s_p[1].c = 2; /* P0 */
244: s_p[2].i = i+1;s_p[2].j = j+1;s_p[2].c = 2; /* P1 */
245: s_p[3].i = i+1;s_p[3].j = j;s_p[3].c = 2; /* P1 */
246: return(0);
247: }
251: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
252: {
254: PetscMPIInt rank;
255: PetscInt proc_I,proc_J;
256: PetscInt cpu_x,cpu_y;
257: PetscInt local_mx,local_my;
258: Vec vlx,vly;
259: PetscInt *LX,*LY,i;
260: PetscScalar *_a;
261: Vec V_SEQ;
262: VecScatter ctx;
265: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
267: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
269: proc_J = rank/cpu_x;
270: proc_I = rank-cpu_x*proc_J;
272: PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);
273: PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);
275: DMDAGetLocalElementSize(da,&local_mx,&local_my,NULL);
276: VecCreate(PETSC_COMM_WORLD,&vlx);
277: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
278: VecSetFromOptions(vlx);
280: VecCreate(PETSC_COMM_WORLD,&vly);
281: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
282: VecSetFromOptions(vly);
284: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
285: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
286: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
287: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
291: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
292: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
293: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
294: VecGetArray(V_SEQ,&_a);
295: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
296: VecRestoreArray(V_SEQ,&_a);
297: VecScatterDestroy(&ctx);
298: VecDestroy(&V_SEQ);
300: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
301: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
302: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
303: VecGetArray(V_SEQ,&_a);
304: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
305: VecRestoreArray(V_SEQ,&_a);
306: VecScatterDestroy(&ctx);
307: VecDestroy(&V_SEQ);
311: *_lx = LX;
312: *_ly = LY;
314: VecDestroy(&vlx);
315: VecDestroy(&vly);
316: return(0);
317: }
321: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
322: {
323: DM cda;
324: Vec coords;
325: DMDACoor2d **_coords;
326: PetscInt si,sj,nx,ny,i,j;
327: FILE *fp;
328: char fname[PETSC_MAX_PATH_LEN];
329: PetscMPIInt rank;
333: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
334: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
335: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
336: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
338: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
340: DMGetCoordinateDM(da,&cda);
341: DMGetCoordinatesLocal(da,&coords);
342: DMDAVecGetArray(cda,coords,&_coords);
343: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
344: for (j = sj; j < sj+ny-1; j++) {
345: for (i = si; i < si+nx-1; i++) {
346: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
347: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i].x),(double)PetscRealPart(_coords[j+1][i].y));
348: 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));
349: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i+1].x),(double)PetscRealPart(_coords[j][i+1].y));
350: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
351: }
352: }
353: DMDAVecRestoreArray(cda,coords,&_coords);
355: PetscFClose(PETSC_COMM_SELF,fp);
356: return(0);
357: }
361: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
362: {
363: DM cda;
364: Vec coords,local_fields;
365: DMDACoor2d **_coords;
366: FILE *fp;
367: char fname[PETSC_MAX_PATH_LEN];
368: PetscMPIInt rank;
369: PetscInt si,sj,nx,ny,i,j;
370: PetscInt n_dofs,d;
371: PetscScalar *_fields;
375: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
376: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
377: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
378: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
380: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
381: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
382: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
383: for (d = 0; d < n_dofs; d++) {
384: const char *field_name;
385: DMDAGetFieldName(da,d,&field_name);
386: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
387: }
388: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
391: DMGetCoordinateDM(da,&cda);
392: DMGetCoordinatesLocal(da,&coords);
393: DMDAVecGetArray(cda,coords,&_coords);
394: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
396: DMCreateLocalVector(da,&local_fields);
397: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
398: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
399: VecGetArray(local_fields,&_fields);
402: for (j = sj; j < sj+ny; j++) {
403: for (i = si; i < si+nx; i++) {
404: PetscScalar coord_x,coord_y;
405: PetscScalar field_d;
407: coord_x = _coords[j][i].x;
408: coord_y = _coords[j][i].y;
410: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
411: for (d = 0; d < n_dofs; d++) {
412: field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
413: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
414: }
415: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
416: }
417: }
418: VecRestoreArray(local_fields,&_fields);
419: VecDestroy(&local_fields);
421: DMDAVecRestoreArray(cda,coords,&_coords);
423: PetscFClose(PETSC_COMM_SELF,fp);
424: return(0);
425: }
429: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
430: {
431: DM cda;
432: Vec local_fields;
433: FILE *fp;
434: char fname[PETSC_MAX_PATH_LEN];
435: PetscMPIInt rank;
436: PetscInt si,sj,nx,ny,i,j,p;
437: PetscInt n_dofs,d;
438: GaussPointCoefficients **_coefficients;
439: PetscErrorCode ierr;
442: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
443: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
444: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
445: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
447: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
448: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
449: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
450: for (d = 0; d < n_dofs; d++) {
451: const char *field_name;
452: DMDAGetFieldName(da,d,&field_name);
453: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
454: }
455: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
458: DMGetCoordinateDM(da,&cda);
459: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
461: DMCreateLocalVector(da,&local_fields);
462: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
463: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
464: DMDAVecGetArray(da,local_fields,&_coefficients);
467: for (j = sj; j < sj+ny; j++) {
468: for (i = si; i < si+nx; i++) {
469: PetscScalar coord_x,coord_y;
471: for (p = 0; p < GAUSS_POINTS; p++) {
472: coord_x = _coefficients[j][i].gp_coords[2*p];
473: coord_y = _coefficients[j][i].gp_coords[2*p+1];
475: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
477: 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]));
478: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
479: }
480: }
481: }
482: DMDAVecRestoreArray(da,local_fields,&_coefficients);
483: VecDestroy(&local_fields);
485: PetscFClose(PETSC_COMM_SELF,fp);
486: return(0);
487: }
490: 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)
491: {
492: PetscInt ij;
493: PetscInt r,c,nc;
495: nc = u_NPE*u_dof;
497: r = w_dof*wi+wd;
498: c = u_dof*ui+ud;
500: ij = r*nc+c;
502: return ij;
503: }
505: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
506: {
507: PetscInt ngp;
508: PetscScalar gp_xi[GAUSS_POINTS][2];
509: PetscScalar gp_weight[GAUSS_POINTS];
510: PetscInt p,i,j,k;
511: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
512: PetscScalar J_p,tildeD[3];
513: PetscScalar B[3][U_DOFS*NODES_PER_EL];
516: /* define quadrature rule */
517: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
519: /* evaluate integral */
520: for (p = 0; p < ngp; p++) {
521: ConstructQ12D_GNi(gp_xi[p],GNi_p);
522: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
524: for (i = 0; i < NODES_PER_EL; i++) {
525: PetscScalar d_dx_i = GNx_p[0][i];
526: PetscScalar d_dy_i = GNx_p[1][i];
528: B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
529: B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
530: B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
531: }
534: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
535: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
536: tildeD[2] = gp_weight[p]*J_p*eta[p];
538: /* form Bt tildeD B */
539: /*
540: Ke_ij = Bt_ik . D_kl . B_lj
541: = B_ki . D_kl . B_lj
542: = B_ki . D_kk . B_kj
543: */
544: for (i = 0; i < 8; i++) {
545: for (j = 0; j < 8; j++) {
546: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
547: Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
548: }
549: }
550: }
551: }
552: }
554: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
555: {
556: PetscInt ngp;
557: PetscScalar gp_xi[GAUSS_POINTS][2];
558: PetscScalar gp_weight[GAUSS_POINTS];
559: PetscInt p,i,j,di;
560: PetscScalar Ni_p[NODES_PER_EL];
561: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
562: PetscScalar J_p,fac;
565: /* define quadrature rule */
566: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
568: /* evaluate integral */
569: for (p = 0; p < ngp; p++) {
570: ConstructQ12D_Ni(gp_xi[p],Ni_p);
571: ConstructQ12D_GNi(gp_xi[p],GNi_p);
572: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
573: fac = gp_weight[p]*J_p;
575: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
576: for (di = 0; di < NSD; di++) { /* u dofs */
577: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
578: PetscInt IJ;
579: /* Ke[4*u_idx+j] = Ke[4*u_idx+j] - GNx_p[di][i] * Ni_p[j] * fac; */
580: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);
582: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
583: }
584: }
585: }
586: }
587: }
589: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
590: {
591: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
592: PetscInt i,j;
593: PetscInt nr_g,nc_g;
595: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
596: FormGradientOperatorQ1(Ge,coords);
598: nr_g = U_DOFS*NODES_PER_EL;
599: nc_g = P_DOFS*NODES_PER_EL;
601: for (i = 0; i < nr_g; i++) {
602: for (j = 0; j < nc_g; j++) {
603: De[nr_g*j+i] = Ge[nc_g*i+j];
604: }
605: }
606: }
608: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
609: {
610: PetscInt ngp;
611: PetscScalar gp_xi[GAUSS_POINTS][2];
612: PetscScalar gp_weight[GAUSS_POINTS];
613: PetscInt p,i,j;
614: PetscScalar Ni_p[NODES_PER_EL];
615: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
616: PetscScalar J_p,fac,eta_avg;
619: /* define quadrature rule */
620: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
622: /* evaluate integral */
623: for (p = 0; p < ngp; p++) {
624: ConstructQ12D_Ni(gp_xi[p],Ni_p);
625: ConstructQ12D_GNi(gp_xi[p],GNi_p);
626: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
627: fac = gp_weight[p]*J_p;
629: for (i = 0; i < NODES_PER_EL; i++) {
630: for (j = 0; j < NODES_PER_EL; j++) {
631: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
632: }
633: }
634: }
636: /* scale */
637: eta_avg = 0.0;
638: for (p = 0; p < ngp; p++) eta_avg += eta[p];
639: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
640: fac = 1.0/eta_avg;
641: for (i = 0; i < NODES_PER_EL; i++) {
642: for (j = 0; j < NODES_PER_EL; j++) {
643: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
644: }
645: }
646: }
648: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
649: {
650: PetscInt ngp;
651: PetscScalar gp_xi[GAUSS_POINTS][2];
652: PetscScalar gp_weight[GAUSS_POINTS];
653: PetscInt p,i,j;
654: PetscScalar Ni_p[NODES_PER_EL];
655: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
656: PetscScalar J_p,fac,eta_avg;
659: /* define quadrature rule */
660: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
662: /* evaluate integral */
663: for (p = 0; p < ngp; p++) {
664: ConstructQ12D_Ni(gp_xi[p],Ni_p);
665: ConstructQ12D_GNi(gp_xi[p],GNi_p);
666: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
667: fac = gp_weight[p]*J_p;
669: for (i = 0; i < NODES_PER_EL; i++) {
670: for (j = 0; j < NODES_PER_EL; j++) {
671: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
672: }
673: }
674: }
676: /* scale */
677: eta_avg = 0.0;
678: for (p = 0; p < ngp; p++) eta_avg += eta[p];
679: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
680: fac = 1.0/eta_avg;
681: for (i = 0; i < NODES_PER_EL; i++) {
682: for (j = 0; j < NODES_PER_EL; j++) {
683: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
684: }
685: }
686: }
688: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
689: {
690: PetscInt ngp;
691: PetscScalar gp_xi[GAUSS_POINTS][2];
692: PetscScalar gp_weight[GAUSS_POINTS];
693: PetscInt p,i;
694: PetscScalar Ni_p[NODES_PER_EL];
695: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
696: PetscScalar J_p,fac;
699: /* define quadrature rule */
700: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
702: /* evaluate integral */
703: for (p = 0; p < ngp; p++) {
704: ConstructQ12D_Ni(gp_xi[p],Ni_p);
705: ConstructQ12D_GNi(gp_xi[p],GNi_p);
706: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
707: fac = gp_weight[p]*J_p;
709: for (i = 0; i < NODES_PER_EL; i++) {
710: Fe[NSD*i] += fac*Ni_p[i]*fx[p];
711: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
712: }
713: }
714: }
718: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
719: {
721: /* get coords for the element */
722: el_coords[NSD*0+0] = _coords[ej][ei].x;el_coords[NSD*0+1] = _coords[ej][ei].y;
723: el_coords[NSD*1+0] = _coords[ej+1][ei].x;el_coords[NSD*1+1] = _coords[ej+1][ei].y;
724: el_coords[NSD*2+0] = _coords[ej+1][ei+1].x;el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
725: el_coords[NSD*3+0] = _coords[ej][ei+1].x;el_coords[NSD*3+1] = _coords[ej][ei+1].y;
726: return(0);
727: }
731: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
732: {
733: DM cda;
734: Vec coords;
735: DMDACoor2d **_coords;
736: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
737: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
738: PetscInt sex,sey,mx,my;
739: PetscInt ei,ej;
740: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
741: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
742: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
743: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
744: PetscScalar el_coords[NODES_PER_EL*NSD];
745: Vec local_properties;
746: GaussPointCoefficients **props;
747: PetscScalar *prop_eta;
748: PetscErrorCode ierr;
751: /* setup for coords */
752: DMGetCoordinateDM(stokes_da,&cda);
753: DMGetCoordinatesLocal(stokes_da,&coords);
754: DMDAVecGetArray(cda,coords,&_coords);
756: /* setup for coefficients */
757: DMCreateLocalVector(properties_da,&local_properties);
758: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
759: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
760: DMDAVecGetArray(properties_da,local_properties,&props);
762: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
763: for (ej = sey; ej < sey+my; ej++) {
764: for (ei = sex; ei < sex+mx; ei++) {
765: /* get coords for the element */
766: GetElementCoords(_coords,ei,ej,el_coords);
768: /* get coefficients for the element */
769: prop_eta = props[ej][ei].eta;
771: /* initialise element stiffness matrix */
772: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
773: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
774: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
775: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
777: /* form element stiffness matrix */
778: FormStressOperatorQ1(Ae,el_coords,prop_eta);
779: FormGradientOperatorQ1(Ge,el_coords);
780: FormDivergenceOperatorQ1(De,el_coords);
781: FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);
783: /* insert element matrix into global matrix */
784: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
785: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
786: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
787: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
788: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
789: }
790: }
791: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
792: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
794: DMDAVecRestoreArray(cda,coords,&_coords);
796: DMDAVecRestoreArray(properties_da,local_properties,&props);
797: VecDestroy(&local_properties);
798: return(0);
799: }
803: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
804: {
805: DM cda;
806: Vec coords;
807: DMDACoor2d **_coords;
808: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
809: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
810: PetscInt sex,sey,mx,my;
811: PetscInt ei,ej;
812: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
813: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
814: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
815: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
816: PetscScalar el_coords[NODES_PER_EL*NSD];
817: Vec local_properties;
818: GaussPointCoefficients **props;
819: PetscScalar *prop_eta;
820: PetscErrorCode ierr;
823: /* setup for coords */
824: DMGetCoordinateDM(stokes_da,&cda);
825: DMGetCoordinatesLocal(stokes_da,&coords);
826: DMDAVecGetArray(cda,coords,&_coords);
828: /* setup for coefficients */
829: DMCreateLocalVector(properties_da,&local_properties);
830: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
831: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
832: DMDAVecGetArray(properties_da,local_properties,&props);
834: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
835: for (ej = sey; ej < sey+my; ej++) {
836: for (ei = sex; ei < sex+mx; ei++) {
837: /* get coords for the element */
838: GetElementCoords(_coords,ei,ej,el_coords);
840: /* get coefficients for the element */
841: prop_eta = props[ej][ei].eta;
843: /* initialise element stiffness matrix */
844: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
845: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
846: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
847: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
850: /* form element stiffness matrix */
851: FormStressOperatorQ1(Ae,el_coords,prop_eta);
852: FormGradientOperatorQ1(Ge,el_coords);
853: /* FormDivergenceOperatorQ1(De, el_coords); */
854: FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);
856: /* insert element matrix into global matrix */
857: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
858: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
859: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
860: /* MatSetValuesStencil(A, NODES_PER_EL*P_DOFS,p_eqn, NODES_PER_EL*U_DOFS,u_eqn, De, ADD_VALUES); */
861: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
862: }
863: }
864: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
865: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
867: DMDAVecRestoreArray(cda,coords,&_coords);
869: DMDAVecRestoreArray(properties_da,local_properties,&props);
870: VecDestroy(&local_properties);
871: return(0);
872: }
876: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
877: {
878: PetscInt n;
881: for (n = 0; n < 4; n++) {
882: 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];
883: 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];
884: 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];
885: }
886: return(0);
887: }
891: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
892: {
893: DM cda;
894: Vec coords;
895: DMDACoor2d **_coords;
896: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
897: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
898: PetscInt sex,sey,mx,my;
899: PetscInt ei,ej;
900: PetscScalar Fe[NODES_PER_EL*U_DOFS];
901: PetscScalar He[NODES_PER_EL*P_DOFS];
902: PetscScalar el_coords[NODES_PER_EL*NSD];
903: Vec local_properties;
904: GaussPointCoefficients **props;
905: PetscScalar *prop_fx,*prop_fy;
906: Vec local_F;
907: StokesDOF **ff;
908: PetscErrorCode ierr;
911: /* setup for coords */
912: DMGetCoordinateDM(stokes_da,&cda);
913: DMGetCoordinatesLocal(stokes_da,&coords);
914: DMDAVecGetArray(cda,coords,&_coords);
916: /* setup for coefficients */
917: DMGetLocalVector(properties_da,&local_properties);
918: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
919: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
920: DMDAVecGetArray(properties_da,local_properties,&props);
922: /* get acces to the vector */
923: DMGetLocalVector(stokes_da,&local_F);
924: VecZeroEntries(local_F);
925: DMDAVecGetArray(stokes_da,local_F,&ff);
928: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
929: for (ej = sey; ej < sey+my; ej++) {
930: for (ei = sex; ei < sex+mx; ei++) {
931: /* get coords for the element */
932: GetElementCoords(_coords,ei,ej,el_coords);
934: /* get coefficients for the element */
935: prop_fx = props[ej][ei].fx;
936: prop_fy = props[ej][ei].fy;
938: /* initialise element stiffness matrix */
939: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
940: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
943: /* form element stiffness matrix */
944: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
946: /* insert element matrix into global matrix */
947: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
949: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
950: }
951: }
953: DMDAVecRestoreArray(stokes_da,local_F,&ff);
954: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
955: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
956: DMRestoreLocalVector(stokes_da,&local_F);
959: DMDAVecRestoreArray(cda,coords,&_coords);
961: DMDAVecRestoreArray(properties_da,local_properties,&props);
962: DMRestoreLocalVector(properties_da,&local_properties);
963: return(0);
964: }
968: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM *_da,Vec *_X)
969: {
970: DM da,cda;
971: Vec X,local_X;
972: StokesDOF **_stokes;
973: Vec coords;
974: DMDACoor2d **_coords;
975: PetscInt si,sj,ei,ej,i,j;
979: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
980: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,&da);
981: DMDASetFieldName(da,0,"anlytic_Vx");
982: DMDASetFieldName(da,1,"anlytic_Vy");
983: DMDASetFieldName(da,2,"analytic_P");
986: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.,0.);
989: DMGetCoordinatesLocal(da,&coords);
990: DMGetCoordinateDM(da,&cda);
991: DMDAVecGetArray(cda,coords,&_coords);
993: DMCreateGlobalVector(da,&X);
994: DMCreateLocalVector(da,&local_X);
995: DMDAVecGetArray(da,local_X,&_stokes);
997: DMDAGetGhostCorners(da,&si,&sj,0,&ei,&ej,0);
998: for (j = sj; j < sj+ej; j++) {
999: for (i = si; i < si+ei; i++) {
1000: PetscReal pos[2],pressure,vel[2],total_stress[3],strain_rate[3];
1002: pos[0] = PetscRealPart(_coords[j][i].x);
1003: pos[1] = PetscRealPart(_coords[j][i].y);
1005: evaluate_solCx(pos,eta0,eta1,xc,nz,vel,&pressure,total_stress,strain_rate);
1007: _stokes[j][i].u_dof = vel[0];
1008: _stokes[j][i].v_dof = vel[1];
1009: _stokes[j][i].p_dof = pressure;
1010: }
1011: }
1012: DMDAVecRestoreArray(da,local_X,&_stokes);
1013: DMDAVecRestoreArray(cda,coords,&_coords);
1015: DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1016: DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);
1018: VecDestroy(&local_X);
1020: *_da = da;
1021: *_X = X;
1022: return(0);
1023: }
1027: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
1028: {
1030: /* get the nodal fields */
1031: 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;
1032: 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;
1033: 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;
1034: 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;
1035: return(0);
1036: }
1040: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
1041: {
1042: DM cda;
1043: Vec coords,X_analytic_local,X_local;
1044: DMDACoor2d **_coords;
1045: PetscInt sex,sey,mx,my;
1046: PetscInt ei,ej;
1047: PetscScalar el_coords[NODES_PER_EL*NSD];
1048: StokesDOF **stokes_analytic,**stokes;
1049: StokesDOF stokes_analytic_e[4],stokes_e[4];
1051: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1052: PetscScalar Ni_p[NODES_PER_EL];
1053: PetscInt ngp;
1054: PetscScalar gp_xi[GAUSS_POINTS][2];
1055: PetscScalar gp_weight[GAUSS_POINTS];
1056: PetscInt p,i;
1057: PetscScalar J_p,fac;
1058: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1059: PetscInt M;
1060: PetscReal xymin[2],xymax[2];
1064: /* define quadrature rule */
1065: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1067: /* setup for coords */
1068: DMGetCoordinateDM(stokes_da,&cda);
1069: DMGetCoordinatesLocal(stokes_da,&coords);
1070: DMDAVecGetArray(cda,coords,&_coords);
1072: /* setup for analytic */
1073: DMCreateLocalVector(stokes_da,&X_analytic_local);
1074: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1075: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1076: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1078: /* setup for solution */
1079: DMCreateLocalVector(stokes_da,&X_local);
1080: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1081: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1082: DMDAVecGetArray(stokes_da,X_local,&stokes);
1084: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1085: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1087: h = (xymax[0]-xymin[0])/((PetscReal)M);
1089: tp_L2 = tu_L2 = tu_H1 = 0.0;
1091: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
1092: for (ej = sey; ej < sey+my; ej++) {
1093: for (ei = sex; ei < sex+mx; ei++) {
1094: /* get coords for the element */
1095: GetElementCoords(_coords,ei,ej,el_coords);
1096: StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1097: StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);
1099: /* evaluate integral */
1100: p_e_L2 = 0.0;
1101: u_e_L2 = 0.0;
1102: u_e_H1 = 0.0;
1103: for (p = 0; p < ngp; p++) {
1104: ConstructQ12D_Ni(gp_xi[p],Ni_p);
1105: ConstructQ12D_GNi(gp_xi[p],GNi_p);
1106: ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1107: fac = gp_weight[p]*J_p;
1109: for (i = 0; i < NODES_PER_EL; i++) {
1110: PetscScalar u_error,v_error;
1112: 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);
1114: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1115: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1116: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);
1118: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1119: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error /* du/dy */
1120: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1121: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error); /* dv/dy */
1122: }
1123: }
1125: tp_L2 += p_e_L2;
1126: tu_L2 += u_e_L2;
1127: tu_H1 += u_e_H1;
1128: }
1129: }
1130: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1131: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1132: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1133: p_L2 = PetscSqrtScalar(p_L2);
1134: u_L2 = PetscSqrtScalar(u_L2);
1135: u_H1 = PetscSqrtScalar(u_H1);
1137: 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));
1140: DMDAVecRestoreArray(cda,coords,&_coords);
1142: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1143: VecDestroy(&X_analytic_local);
1144: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1145: VecDestroy(&X_local);
1146: return(0);
1147: }
1151: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1152: {
1153: DM da_Stokes,da_prop;
1154: PetscInt u_dof,p_dof,dof,stencil_width;
1155: Mat A,B;
1156: PetscInt mxl,myl;
1157: DM prop_cda,vel_cda;
1158: Vec prop_coords,vel_coords;
1159: PetscInt si,sj,nx,ny,i,j,p;
1160: Vec f,X;
1161: PetscInt prop_dof,prop_stencil_width;
1162: Vec properties,l_properties;
1163: PetscReal dx,dy;
1164: PetscInt M,N;
1165: DMDACoor2d **_prop_coords,**_vel_coords;
1166: GaussPointCoefficients **element_props;
1167: PetscInt its;
1168: KSP ksp_S;
1169: PetscInt coefficient_structure = 0;
1170: PetscInt cpu_x,cpu_y,*lx = NULL,*ly = NULL;
1171: PetscBool use_gp_coords = PETSC_FALSE,set;
1172: char filename[PETSC_MAX_PATH_LEN];
1173: PetscErrorCode ierr;
1176: /* Generate the da for velocity and pressure */
1177: /*
1178: We use Q1 elements for the temperature.
1179: FEM has a 9-point stencil (BOX) or connectivity pattern
1180: Num nodes in each direction is mx+1, my+1
1181: */
1182: u_dof = U_DOFS; /* Vx, Vy - velocities */
1183: p_dof = P_DOFS; /* p - pressure */
1184: dof = u_dof+p_dof;
1185: stencil_width = 1;
1186: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1187: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&da_Stokes);
1188: DMDASetFieldName(da_Stokes,0,"Vx");
1189: DMDASetFieldName(da_Stokes,1,"Vy");
1190: DMDASetFieldName(da_Stokes,2,"P");
1192: /* unit box [0,1] x [0,1] */
1193: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.,0.);
1196: /* Generate element properties, we will assume all material properties are constant over the element */
1197: /* local number of elements */
1198: DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,NULL);
1200: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1201: DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1202: DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);
1204: prop_dof = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1205: prop_stencil_width = 0;
1206: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1207: mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
1208: PetscFree(lx);
1209: PetscFree(ly);
1211: /* define centroid positions */
1212: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1213: dx = 1.0/((PetscReal)(M));
1214: dy = 1.0/((PetscReal)(N));
1216: 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);
1218: /* define coefficients */
1219: PetscOptionsGetInt(NULL,"-c_str",&coefficient_structure,NULL);
1220: /* PetscPrintf(PETSC_COMM_WORLD, "Using coeficient structure %D \n", coefficient_structure); */
1222: DMCreateGlobalVector(da_prop,&properties);
1223: DMCreateLocalVector(da_prop,&l_properties);
1224: DMDAVecGetArray(da_prop,l_properties,&element_props);
1226: DMGetCoordinateDM(da_prop,&prop_cda);
1227: DMGetCoordinatesLocal(da_prop,&prop_coords);
1228: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
1230: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
1232: DMGetCoordinateDM(da_Stokes,&vel_cda);
1233: DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1234: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1237: /* interpolate the coordinates */
1238: for (j = sj; j < sj+ny; j++) {
1239: for (i = si; i < si+nx; i++) {
1240: PetscInt ngp;
1241: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1242: PetscScalar el_coords[8];
1244: GetElementCoords(_vel_coords,i,j,el_coords);
1245: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1247: for (p = 0; p < GAUSS_POINTS; p++) {
1248: PetscScalar gp_x,gp_y;
1249: PetscInt n;
1250: PetscScalar xi_p[2],Ni_p[4];
1252: xi_p[0] = gp_xi[p][0];
1253: xi_p[1] = gp_xi[p][1];
1254: ConstructQ12D_Ni(xi_p,Ni_p);
1256: gp_x = 0.0;
1257: gp_y = 0.0;
1258: for (n = 0; n < NODES_PER_EL; n++) {
1259: gp_x = gp_x+Ni_p[n]*el_coords[2*n];
1260: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1261: }
1262: element_props[j][i].gp_coords[2*p] = gp_x;
1263: element_props[j][i].gp_coords[2*p+1] = gp_y;
1264: }
1265: }
1266: }
1268: /* define the coefficients */
1269: PetscOptionsGetBool(NULL,"-use_gp_coords",&use_gp_coords,0);
1271: for (j = sj; j < sj+ny; j++) {
1272: for (i = si; i < si+nx; i++) {
1273: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1274: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1275: PetscReal coord_x,coord_y;
1277: if (coefficient_structure == 0) {
1278: PetscReal opts_eta0,opts_eta1,opts_xc;
1279: PetscInt opts_nz;
1281: opts_eta0 = 1.0;
1282: opts_eta1 = 1.0;
1283: opts_xc = 0.5;
1284: opts_nz = 1;
1286: PetscOptionsGetReal(NULL,"-solcx_eta0",&opts_eta0,0);
1287: PetscOptionsGetReal(NULL,"-solcx_eta1",&opts_eta1,0);
1288: PetscOptionsGetReal(NULL,"-solcx_xc",&opts_xc,0);
1289: PetscOptionsGetInt(NULL,"-solcx_nz",&opts_nz,0);
1291: for (p = 0; p < GAUSS_POINTS; p++) {
1292: coord_x = centroid_x;
1293: coord_y = centroid_y;
1294: if (use_gp_coords) {
1295: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1296: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1297: }
1300: element_props[j][i].eta[p] = opts_eta0;
1301: if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;
1303: element_props[j][i].fx[p] = 0.0;
1304: element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1305: }
1306: } else if (coefficient_structure == 1) { /* square sinker */
1307: PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;
1309: opts_eta0 = 1.0;
1310: opts_eta1 = 1.0;
1311: opts_dx = 0.50;
1312: opts_dy = 0.50;
1314: PetscOptionsGetReal(NULL,"-sinker_eta0",&opts_eta0,0);
1315: PetscOptionsGetReal(NULL,"-sinker_eta1",&opts_eta1,0);
1316: PetscOptionsGetReal(NULL,"-sinker_dx",&opts_dx,0);
1317: PetscOptionsGetReal(NULL,"-sinker_dy",&opts_dy,0);
1320: for (p = 0; p < GAUSS_POINTS; p++) {
1321: coord_x = centroid_x;
1322: coord_y = centroid_y;
1323: if (use_gp_coords) {
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: }
1328: element_props[j][i].eta[p] = opts_eta0;
1329: element_props[j][i].fx[p] = 0.0;
1330: element_props[j][i].fy[p] = 0.0;
1332: if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1333: if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1334: element_props[j][i].eta[p] = opts_eta1;
1335: element_props[j][i].fx[p] = 0.0;
1336: element_props[j][i].fy[p] = -1.0;
1337: }
1338: }
1339: }
1340: } else if (coefficient_structure == 2) { /* circular sinker */
1341: PetscReal opts_eta0,opts_eta1,opts_r,radius2;
1343: opts_eta0 = 1.0;
1344: opts_eta1 = 1.0;
1345: opts_r = 0.25;
1347: PetscOptionsGetReal(NULL,"-sinker_eta0",&opts_eta0,0);
1348: PetscOptionsGetReal(NULL,"-sinker_eta1",&opts_eta1,0);
1349: PetscOptionsGetReal(NULL,"-sinker_r",&opts_r,0);
1351: for (p = 0; p < GAUSS_POINTS; p++) {
1352: coord_x = centroid_x;
1353: coord_y = centroid_y;
1354: if (use_gp_coords) {
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: radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1364: if (radius2 < opts_r*opts_r) {
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: } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1371: PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;
1373: opts_eta0 = 1.0;
1374: opts_eta1 = 1.0;
1375: opts_r = 0.25;
1376: opts_c0x = 0.35; /* circle center */
1377: opts_c0y = 0.35;
1378: opts_s0x = 0.7; /* square center */
1379: opts_s0y = 0.7;
1380: opts_dx = 0.25;
1381: opts_dy = 0.25;
1382: opts_phi = 25;
1384: PetscOptionsGetReal(NULL,"-sinker_eta0",&opts_eta0,0);
1385: PetscOptionsGetReal(NULL,"-sinker_eta1",&opts_eta1,0);
1386: PetscOptionsGetReal(NULL,"-sinker_r",&opts_r,0);
1387: PetscOptionsGetReal(NULL,"-sinker_c0x",&opts_c0x,0);
1388: PetscOptionsGetReal(NULL,"-sinker_c0y",&opts_c0y,0);
1389: PetscOptionsGetReal(NULL,"-sinker_s0x",&opts_s0x,0);
1390: PetscOptionsGetReal(NULL,"-sinker_s0y",&opts_s0y,0);
1391: PetscOptionsGetReal(NULL,"-sinker_dx",&opts_dx,0);
1392: PetscOptionsGetReal(NULL,"-sinker_dy",&opts_dy,0);
1393: PetscOptionsGetReal(NULL,"-sinker_phi",&opts_phi,0);
1394: opts_phi *= PETSC_PI / 180;
1396: for (p = 0; p < GAUSS_POINTS; p++) {
1397: coord_x = centroid_x;
1398: coord_y = centroid_y;
1399: if (use_gp_coords) {
1400: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1401: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1402: }
1404: element_props[j][i].eta[p] = opts_eta0;
1405: element_props[j][i].fx[p] = 0.0;
1406: element_props[j][i].fy[p] = -0.2;
1408: radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1409: if (radius2 < opts_r*opts_r
1410: || (PetscAbs(+(coord_x - opts_s0x)*PetscCosReal(opts_phi) + (coord_y - opts_s0y)*PetscSinReal(opts_phi)) < opts_dx/2
1411: && PetscAbs(-(coord_x - opts_s0x)*PetscSinReal(opts_phi) + (coord_y - opts_s0y)*PetscCosReal(opts_phi)) < opts_dy/2)) {
1412: element_props[j][i].eta[p] = opts_eta1;
1413: element_props[j][i].fx[p] = 0.0;
1414: element_props[j][i].fy[p] = -1.0;
1415: }
1416: }
1417: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1418: }
1419: }
1420: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
1422: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1424: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1425: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1426: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
1429: DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1430: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for Stokes eqn.","properties");
1433: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1434: DMSetMatType(da_Stokes,MATAIJ);
1435: DMCreateMatrix(da_Stokes,&A);
1436: DMCreateMatrix(da_Stokes,&B);
1437: DMCreateGlobalVector(da_Stokes,&f);
1438: DMCreateGlobalVector(da_Stokes,&X);
1440: /* assemble A11 */
1441: MatZeroEntries(A);
1442: MatZeroEntries(B);
1443: VecZeroEntries(f);
1445: AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1446: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1447: /* build force vector */
1448: AssembleF_Stokes(f,da_Stokes,da_prop,properties);
1450: DMDABCApplyFreeSlip(da_Stokes,A,f);
1451: DMDABCApplyFreeSlip(da_Stokes,B,NULL);
1453: /* SOLVE */
1454: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1455: KSPSetOptionsPrefix(ksp_S,"stokes_");
1456: KSPSetOperators(ksp_S,A,B);
1457: KSPSetDM(ksp_S,da_Stokes);
1458: KSPSetDMActive(ksp_S,PETSC_FALSE);
1459: KSPSetFromOptions(ksp_S);
1460: {
1461: PC pc;
1462: const PetscInt ufields[] = {0,1},pfields[1] = {2};
1463: KSPGetPC(ksp_S,&pc);
1464: PCFieldSplitSetBlockSize(pc,3);
1465: PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1466: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1467: }
1469: KSPSolve(ksp_S,f,X);
1471: PetscOptionsGetString(NULL,"-o",filename,sizeof(filename),&set);
1472: if (set) {
1473: char *ext;
1474: PetscViewer viewer;
1475: PetscBool flg;
1476: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1477: PetscStrrchr(filename,'.',&ext);
1478: PetscStrcmp("vts",ext,&flg);
1479: if (flg) {
1480: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1481: } else {
1482: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1483: }
1484: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1485: PetscViewerFileSetName(viewer,filename);
1486: VecView(X,viewer);
1487: PetscViewerDestroy(&viewer);
1488: }
1489: DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
1491: KSPGetIterationNumber(ksp_S,&its);
1493: if (coefficient_structure == 0) {
1494: PetscReal opts_eta0,opts_eta1,opts_xc;
1495: PetscInt opts_nz,N;
1496: DM da_Stokes_analytic;
1497: Vec X_analytic;
1498: PetscReal nrm1[3],nrm2[3],nrmI[3];
1500: opts_eta0 = 1.0;
1501: opts_eta1 = 1.0;
1502: opts_xc = 0.5;
1503: opts_nz = 1;
1505: PetscOptionsGetReal(NULL,"-solcx_eta0",&opts_eta0,0);
1506: PetscOptionsGetReal(NULL,"-solcx_eta1",&opts_eta1,0);
1507: PetscOptionsGetReal(NULL,"-solcx_xc",&opts_xc,0);
1508: PetscOptionsGetInt(NULL,"-solcx_nz",&opts_nz,0);
1511: DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1512: DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
1514: DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);
1517: VecAXPY(X_analytic,-1.0,X);
1518: VecGetSize(X_analytic,&N);
1519: N = N/3;
1521: VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1522: VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1523: VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);
1525: VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1526: VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1527: VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);
1529: VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1530: VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1531: VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);
1533: DMDestroy(&da_Stokes_analytic);
1534: VecDestroy(&X_analytic);
1535: }
1538: KSPDestroy(&ksp_S);
1539: VecDestroy(&X);
1540: VecDestroy(&f);
1541: MatDestroy(&A);
1542: MatDestroy(&B);
1544: DMDestroy(&da_Stokes);
1545: DMDestroy(&da_prop);
1547: VecDestroy(&properties);
1548: VecDestroy(&l_properties);
1549: return(0);
1550: }
1554: int main(int argc,char **args)
1555: {
1557: PetscInt mx,my;
1559: PetscInitialize(&argc,&args,(char*)0,help);
1561: mx = my = 10;
1562: PetscOptionsGetInt(NULL,"-mx",&mx,NULL);
1563: PetscOptionsGetInt(NULL,"-my",&my,NULL);
1565: solve_stokes_2d_coupled(mx,my);
1567: PetscFinalize();
1568: return 0;
1569: }
1571: /* -------------------------- helpers for boundary conditions -------------------------------- */
1575: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1576: {
1577: DM cda;
1578: Vec coords;
1579: PetscInt si,sj,nx,ny,i,j;
1580: PetscInt M,N;
1581: DMDACoor2d **_coords;
1582: const PetscInt *g_idx;
1583: PetscInt *bc_global_ids;
1584: PetscScalar *bc_vals;
1585: PetscInt nbcs;
1586: PetscInt n_dofs;
1587: PetscErrorCode ierr;
1588: ISLocalToGlobalMapping ltogm;
1591: /* enforce bc's */
1592: DMGetLocalToGlobalMapping(da,<ogm);
1593: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1595: DMGetCoordinateDM(da,&cda);
1596: DMGetCoordinatesLocal(da,&coords);
1597: DMDAVecGetArray(cda,coords,&_coords);
1598: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1599: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1601: /* --- */
1603: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1604: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1606: /* init the entries to -1 so VecSetValues will ignore them */
1607: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1609: i = nx-1;
1610: for (j = 0; j < ny; j++) {
1611: PetscInt local_id;
1613: local_id = i+j*nx;
1615: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1617: bc_vals[j] = bc_val;
1618: }
1619: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1620: nbcs = 0;
1621: if ((si+nx) == (M)) nbcs = ny;
1623: if (b != NULL) {
1624: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1625: VecAssemblyBegin(b);
1626: VecAssemblyEnd(b);
1627: }
1628: if (A != NULL) {
1629: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1630: }
1633: PetscFree(bc_vals);
1634: PetscFree(bc_global_ids);
1636: DMDAVecRestoreArray(cda,coords,&_coords);
1637: return(0);
1638: }
1642: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1643: {
1644: DM cda;
1645: Vec coords;
1646: PetscInt si,sj,nx,ny,i,j;
1647: PetscInt M,N;
1648: DMDACoor2d **_coords;
1649: const PetscInt *g_idx;
1650: PetscInt *bc_global_ids;
1651: PetscScalar *bc_vals;
1652: PetscInt nbcs;
1653: PetscInt n_dofs;
1654: PetscErrorCode ierr;
1655: ISLocalToGlobalMapping ltogm;
1658: /* enforce bc's */
1659: DMGetLocalToGlobalMapping(da,<ogm);
1660: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1662: DMGetCoordinateDM(da,&cda);
1663: DMGetCoordinatesLocal(da,&coords);
1664: DMDAVecGetArray(cda,coords,&_coords);
1665: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1666: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1668: /* --- */
1670: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1671: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1673: /* init the entries to -1 so VecSetValues will ignore them */
1674: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1676: i = 0;
1677: for (j = 0; j < ny; j++) {
1678: PetscInt local_id;
1680: local_id = i+j*nx;
1682: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1684: bc_vals[j] = bc_val;
1685: }
1686: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1687: nbcs = 0;
1688: if (si == 0) nbcs = ny;
1690: if (b != NULL) {
1691: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1692: VecAssemblyBegin(b);
1693: VecAssemblyEnd(b);
1694: }
1695: if (A != NULL) {
1696: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1697: }
1700: PetscFree(bc_vals);
1701: PetscFree(bc_global_ids);
1703: DMDAVecRestoreArray(cda,coords,&_coords);
1704: return(0);
1705: }
1709: static PetscErrorCode BCApply_NORTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1710: {
1711: DM cda;
1712: Vec coords;
1713: PetscInt si,sj,nx,ny,i,j;
1714: PetscInt M,N;
1715: DMDACoor2d **_coords;
1716: const PetscInt *g_idx;
1717: PetscInt *bc_global_ids;
1718: PetscScalar *bc_vals;
1719: PetscInt nbcs;
1720: PetscInt n_dofs;
1721: PetscErrorCode ierr;
1722: ISLocalToGlobalMapping ltogm;
1725: /* enforce bc's */
1726: DMGetLocalToGlobalMapping(da,<ogm);
1727: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1729: DMGetCoordinateDM(da,&cda);
1730: DMGetCoordinatesLocal(da,&coords);
1731: DMDAVecGetArray(cda,coords,&_coords);
1732: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1733: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1735: /* --- */
1737: PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1738: PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);
1740: /* init the entries to -1 so VecSetValues will ignore them */
1741: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1743: j = ny-1;
1744: for (i = 0; i < nx; i++) {
1745: PetscInt local_id;
1747: local_id = i+j*nx;
1749: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1751: bc_vals[i] = bc_val;
1752: }
1753: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1754: nbcs = 0;
1755: if ((sj+ny) == (N)) nbcs = nx;
1757: if (b != NULL) {
1758: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1759: VecAssemblyBegin(b);
1760: VecAssemblyEnd(b);
1761: }
1762: if (A != NULL) {
1763: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1764: }
1767: PetscFree(bc_vals);
1768: PetscFree(bc_global_ids);
1770: DMDAVecRestoreArray(cda,coords,&_coords);
1771: return(0);
1772: }
1776: static PetscErrorCode BCApply_SOUTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1777: {
1778: DM cda;
1779: Vec coords;
1780: PetscInt si,sj,nx,ny,i,j;
1781: PetscInt M,N;
1782: DMDACoor2d **_coords;
1783: const PetscInt *g_idx;
1784: PetscInt *bc_global_ids;
1785: PetscScalar *bc_vals;
1786: PetscInt nbcs;
1787: PetscInt n_dofs;
1788: PetscErrorCode ierr;
1789: ISLocalToGlobalMapping ltogm;
1792: /* enforce bc's */
1793: DMGetLocalToGlobalMapping(da,<ogm);
1794: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1796: DMGetCoordinateDM(da,&cda);
1797: DMGetCoordinatesLocal(da,&coords);
1798: DMDAVecGetArray(cda,coords,&_coords);
1799: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1800: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1802: /* --- */
1804: PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1805: PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);
1807: /* init the entries to -1 so VecSetValues will ignore them */
1808: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1810: j = 0;
1811: for (i = 0; i < nx; i++) {
1812: PetscInt local_id;
1814: local_id = i+j*nx;
1816: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1818: bc_vals[i] = bc_val;
1819: }
1820: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1821: nbcs = 0;
1822: if (sj == 0) nbcs = nx;
1824: if (b != NULL) {
1825: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1826: VecAssemblyBegin(b);
1827: VecAssemblyEnd(b);
1828: }
1829: if (A != NULL) {
1830: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1831: }
1834: PetscFree(bc_vals);
1835: PetscFree(bc_global_ids);
1837: DMDAVecRestoreArray(cda,coords,&_coords);
1838: return(0);
1839: }
1843: /*
1844: Free slip sides.
1845: */
1848: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1849: {
1853: BCApply_NORTH(da_Stokes,1,0.0,A,f);
1854: BCApply_EAST(da_Stokes,0,0.0,A,f);
1855: BCApply_SOUTH(da_Stokes,1,0.0,A,f);
1856: BCApply_WEST(da_Stokes,0,0.0,A,f);
1857: return(0);
1858: }