Actual source code: ex43.c
petsc-3.3-p7 2013-05-11
1: static char help[] =
2: "Solves the incompressible, variable viscosity stokes equation in 2d on the unit domain \n\
3: using Q1Q1 elements, stabilized with Bochev's polynomial projection method. \n\
4: The models defined utilise free slip boundary conditions on all sides. \n\
5: Options: \n\
6: -mx : number elements in x-direciton \n\
7: -my : number elements in y-direciton \n\
8: -c_str : indicates the structure of the coefficients to use. \n\
9: -c_str 0 => Setup for an analytic solution with a vertical jump in viscosity. This problem is driven by the \n\
10: forcing function f = ( 0, sin(n_z pi y )cos( pi x ). \n\
11: Parameters: \n\
12: -solcx_eta0 : the viscosity to the left of the interface \n\
13: -solcx_eta1 : the viscosity to the right of the interface \n\
14: -solcx_xc : the location of the interface \n\
15: -solcx_nz : the wavenumber in the y direction \n\
16: -c_str 1 => Setup for a rectangular blob located in the origin of the domain. \n\
17: Parameters: \n\
18: -sinker_eta0 : the viscosity of the background fluid \n\
19: -sinker_eta1 : the viscosity of the blob \n\
20: -sinker_dx : the width of the blob \n\
21: -sinker_dy : the width of the blob \n\
22: -c_str 2 => Setup for a circular blob located in the origin of the domain. \n\
23: Parameters: \n\
24: -sinker_eta0 : the viscosity of the background fluid \n\
25: -sinker_eta1 : the viscosity of the blob \n\
26: -sinker_r : radius of the blob \n\
27: -c_str 3 => Circular and rectangular inclusion\n\
28: Parameters as for cases 1 and 2 (above)\n\
29: -use_gp_coords : evaluate the viscosity and the body force at the global coordinates of the quadrature points.\n\
30: By default, eta and the body force are evaulated at the element center and applied as a constant over the entire element.\n";
32: /* Contributed by Dave May */
34: #include <petscksp.h>
35: #include <petscdmda.h>
37: /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
38: #include "ex43-solCx.h"
40: static PetscErrorCode DMDABCApplyFreeSlip(DM,Mat,Vec);
43: #define NSD 2 /* number of spatial dimensions */
44: #define NODES_PER_EL 4 /* nodes per element */
45: #define U_DOFS 2 /* degrees of freedom per velocity node */
46: #define P_DOFS 1 /* degrees of freedom per pressure node */
47: #define GAUSS_POINTS 4
49: /* cell based evaluation */
50: typedef struct {
51: PetscScalar eta,fx,fy;
52: } Coefficients;
54: /* Gauss point based evaluation 8+4+4+4 = 20 */
55: typedef struct {
56: PetscScalar gp_coords[2*GAUSS_POINTS];
57: PetscScalar eta[GAUSS_POINTS];
58: PetscScalar fx[GAUSS_POINTS];
59: PetscScalar fy[GAUSS_POINTS];
60: } GaussPointCoefficients;
62: typedef struct {
63: PetscScalar u_dof;
64: PetscScalar v_dof;
65: PetscScalar p_dof;
66: } StokesDOF;
69: /*
71: D = [ 2.eta 0 0 ]
72: [ 0 2.eta 0 ]
73: [ 0 0 eta ]
75: B = [ d_dx 0 ]
76: [ 0 d_dy ]
77: [ d_dy d_dx ]
79: */
81: /* FEM routines */
82: /*
83: Element: Local basis function ordering
84: 1-----2
85: | |
86: | |
87: 0-----3
88: */
89: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
90: {
91: PetscScalar xi = _xi[0];
92: PetscScalar eta = _xi[1];
94: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
95: Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
96: Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
97: Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
98: }
100: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
101: {
102: PetscScalar xi = _xi[0];
103: PetscScalar eta = _xi[1];
105: GNi[0][0] = -0.25*(1.0-eta);
106: GNi[0][1] = -0.25*(1.0+eta);
107: GNi[0][2] = 0.25*(1.0+eta);
108: GNi[0][3] = 0.25*(1.0-eta);
110: GNi[1][0] = -0.25*(1.0-xi);
111: GNi[1][1] = 0.25*(1.0-xi);
112: GNi[1][2] = 0.25*(1.0+xi);
113: GNi[1][3] = -0.25*(1.0+xi);
114: }
116: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
117: {
118: PetscScalar J00,J01,J10,J11,J;
119: PetscScalar iJ00,iJ01,iJ10,iJ11;
120: PetscInt i;
122: J00 = J01 = J10 = J11 = 0.0;
123: for (i = 0; i < NODES_PER_EL; i++) {
124: PetscScalar cx = coords[ 2*i+0 ];
125: PetscScalar cy = coords[ 2*i+1 ];
127: J00 = J00+GNi[0][i]*cx; /* J_xx = dx/dxi */
128: J01 = J01+GNi[0][i]*cy; /* J_xy = dy/dxi */
129: J10 = J10+GNi[1][i]*cx; /* J_yx = dx/deta */
130: J11 = J11+GNi[1][i]*cy; /* J_yy = dy/deta */
131: }
132: J = (J00*J11)-(J01*J10);
134: iJ00 = J11/J;
135: iJ01 = -J01/J;
136: iJ10 = -J10/J;
137: iJ11 = J00/J;
140: for (i = 0; i < NODES_PER_EL; i++) {
141: GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
142: GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
143: }
145: if (det_J != NULL) {
146: *det_J = J;
147: }
148: }
150: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
151: {
152: *ngp = 4;
153: gp_xi[0][0] = -0.57735026919;gp_xi[0][1] = -0.57735026919;
154: gp_xi[1][0] = -0.57735026919;gp_xi[1][1] = 0.57735026919;
155: gp_xi[2][0] = 0.57735026919;gp_xi[2][1] = 0.57735026919;
156: gp_xi[3][0] = 0.57735026919;gp_xi[3][1] = -0.57735026919;
157: gp_weight[0] = 1.0;
158: gp_weight[1] = 1.0;
159: gp_weight[2] = 1.0;
160: gp_weight[3] = 1.0;
161: }
164: /* procs to the left claim the ghost node as their element */
167: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
168: {
169: PetscInt m,n,p,M,N,P;
170: PetscInt sx,sy,sz;
173: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
174: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
176: if (mxl != PETSC_NULL) {
177: *mxl = m;
178: if ((sx+m) == M) { /* last proc */
179: *mxl = m-1;
180: }
181: }
182: if (myl != PETSC_NULL) {
183: *myl = n;
184: if ((sy+n) == N) { /* last proc */
185: *myl = n-1;
186: }
187: }
188: if (mzl != PETSC_NULL) {
189: *mzl = p;
190: if ((sz+p) == P) { /* last proc */
191: *mzl = p-1;
192: }
193: }
194: return(0);
195: }
199: static PetscErrorCode DMDAGetElementCorners(DM da,
200: PetscInt *sx,PetscInt *sy,PetscInt *sz,
201: PetscInt *mx,PetscInt *my,PetscInt *mz)
202: {
203: PetscInt si,sj,sk;
206: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
208: *sx = si;
209: if (si != 0) {
210: *sx = si+1;
211: }
213: *sy = sj;
214: if (sj != 0) {
215: *sy = sj+1;
216: }
218: if (sk != PETSC_NULL) {
219: *sz = sk;
220: if (sk != 0) {
221: *sz = sk+1;
222: }
223: }
225: DMDAGetLocalElementSize(da,mx,my,mz);
227: return(0);
228: }
230: /*
231: i,j are the element indices
232: The unknown is a vector quantity.
233: The s[].c is used to indicate the degree of freedom.
234: */
237: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
238: {
240: /* velocity */
241: /* node 0 */
242: s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0; /* Vx0 */
243: s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1; /* Vy0 */
245: /* node 1 */
246: s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0; /* Vx1 */
247: s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1; /* Vy1 */
249: /* node 2 */
250: s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0; /* Vx2 */
251: s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1; /* Vy2 */
253: /* node 3 */
254: s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0; /* Vx3 */
255: s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1; /* Vy3 */
258: /* pressure */
259: s_p[0].i = i;s_p[0].j = j;s_p[0].c = 2; /* P0 */
260: s_p[1].i = i;s_p[1].j = j+1;s_p[1].c = 2; /* P0 */
261: s_p[2].i = i+1;s_p[2].j = j+1;s_p[2].c = 2; /* P1 */
262: s_p[3].i = i+1;s_p[3].j = j;s_p[3].c = 2; /* P1 */
263: return(0);
264: }
268: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
269: {
271: PetscMPIInt rank;
272: PetscInt proc_I,proc_J;
273: PetscInt cpu_x,cpu_y;
274: PetscInt local_mx,local_my;
275: Vec vlx,vly;
276: PetscInt *LX,*LY,i;
277: PetscScalar *_a;
278: Vec V_SEQ;
279: VecScatter ctx;
282: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
284: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
286: proc_J = rank/cpu_x;
287: proc_I = rank-cpu_x*proc_J;
289: PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);
290: PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);
292: DMDAGetLocalElementSize(da,&local_mx,&local_my,PETSC_NULL);
293: VecCreate(PETSC_COMM_WORLD,&vlx);
294: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
295: VecSetFromOptions(vlx);
297: VecCreate(PETSC_COMM_WORLD,&vly);
298: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
299: VecSetFromOptions(vly);
301: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
302: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
303: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
304: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
308: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
309: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
310: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
311: VecGetArray(V_SEQ,&_a);
312: for (i = 0; i < cpu_x; i++) {
313: LX[i] = (PetscInt)PetscRealPart(_a[i]);
314: }
315: VecRestoreArray(V_SEQ,&_a);
316: VecScatterDestroy(&ctx);
317: VecDestroy(&V_SEQ);
319: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
320: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
321: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
322: VecGetArray(V_SEQ,&_a);
323: for (i = 0; i < cpu_y; i++) {
324: LY[i] = (PetscInt)PetscRealPart(_a[i]);
325: }
326: VecRestoreArray(V_SEQ,&_a);
327: VecScatterDestroy(&ctx);
328: VecDestroy(&V_SEQ);
332: *_lx = LX;
333: *_ly = LY;
335: VecDestroy(&vlx);
336: VecDestroy(&vly);
338: return(0);
339: }
343: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
344: {
345: DM cda;
346: Vec coords;
347: DMDACoor2d **_coords;
348: PetscInt si,sj,nx,ny,i,j;
349: FILE *fp;
350: char fname[PETSC_MAX_PATH_LEN];
351: PetscMPIInt rank;
355: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
356: PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
357: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
358: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
360: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
362: DMDAGetCoordinateDA(da,&cda);
363: DMDAGetGhostedCoordinates(da,&coords);
364: DMDAVecGetArray(cda,coords,&_coords);
365: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
366: for (j = sj; j < sj+ny-1; j++) {
367: for (i = si; i < si+nx-1; i++) {
368: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
369: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i].x),PetscRealPart(_coords[j+1][i].y));
370: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i+1].x),PetscRealPart(_coords[j+1][i+1].y));
371: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i+1].x),PetscRealPart(_coords[j][i+1].y));
372: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
373: }
374: }
375: DMDAVecRestoreArray(cda,coords,&_coords);
377: PetscFClose(PETSC_COMM_SELF,fp);
378: return(0);
379: }
383: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
384: {
385: DM cda;
386: Vec coords,local_fields;
387: DMDACoor2d **_coords;
388: FILE *fp;
389: char fname[PETSC_MAX_PATH_LEN];
390: PetscMPIInt rank;
391: PetscInt si,sj,nx,ny,i,j;
392: PetscInt n_dofs,d;
393: PetscScalar *_fields;
397: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
398: PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
399: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
400: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
402: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
403: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
404: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
405: for (d = 0; d < n_dofs; d++) {
406: const char *field_name;
407: DMDAGetFieldName(da,d,&field_name);
408: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
409: }
410: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
413: DMDAGetCoordinateDA(da,&cda);
414: DMDAGetGhostedCoordinates(da,&coords);
415: DMDAVecGetArray(cda,coords,&_coords);
416: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
418: DMCreateLocalVector(da,&local_fields);
419: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
420: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
421: VecGetArray(local_fields,&_fields);
424: for (j = sj; j < sj+ny; j++) {
425: for (i = si; i < si+nx; i++) {
426: PetscScalar coord_x,coord_y;
427: PetscScalar field_d;
429: coord_x = _coords[j][i].x;
430: coord_y = _coords[j][i].y;
432: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
433: for (d = 0; d < n_dofs; d++) {
434: field_d = _fields[ n_dofs*((i-si)+(j-sj)*(nx))+d ];
435: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",PetscRealPart(field_d));
436: }
437: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
438: }
439: }
440: VecRestoreArray(local_fields,&_fields);
441: VecDestroy(&local_fields);
443: DMDAVecRestoreArray(cda,coords,&_coords);
445: PetscFClose(PETSC_COMM_SELF,fp);
446: return(0);
447: }
451: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
452: {
453: DM cda;
454: Vec local_fields;
455: FILE *fp;
456: char fname[PETSC_MAX_PATH_LEN];
457: PetscMPIInt rank;
458: PetscInt si,sj,nx,ny,i,j,p;
459: PetscInt n_dofs,d;
460: GaussPointCoefficients **_coefficients;
461: PetscErrorCode ierr;
464: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
465: PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
466: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
467: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
469: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
470: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
471: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
472: for (d = 0; d < n_dofs; d++) {
473: const char *field_name;
474: DMDAGetFieldName(da,d,&field_name);
475: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
476: }
477: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
480: DMDAGetCoordinateDA(da,&cda);
481: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
483: DMCreateLocalVector(da,&local_fields);
484: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
485: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
486: DMDAVecGetArray(da,local_fields,&_coefficients);
489: for (j = sj; j < sj+ny; j++) {
490: for (i = si; i < si+nx; i++) {
491: PetscScalar coord_x,coord_y;
493: for (p = 0; p < GAUSS_POINTS; p++) {
494: coord_x = _coefficients[j][i].gp_coords[2*p];
495: coord_y = _coefficients[j][i].gp_coords[2*p+1];
497: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
499: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e",PetscRealPart(_coefficients[j][i].eta[p]),PetscRealPart(_coefficients[j][i].fx[p]),PetscRealPart(_coefficients[j][i].fy[p]));
500: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
501: }
502: }
503: }
504: DMDAVecRestoreArray(da,local_fields,&_coefficients);
505: VecDestroy(&local_fields);
507: PetscFClose(PETSC_COMM_SELF,fp);
508: return(0);
509: }
512: static PetscInt ASS_MAP_wIwDI_uJuDJ(
513: PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,
514: PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
515: {
516: PetscInt ij;
517: PetscInt r,c,nc;
519: nc = u_NPE*u_dof;
521: r = w_dof*wi+wd;
522: c = u_dof*ui+ud;
524: ij = r*nc+c;
526: return ij;
527: }
529: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
530: {
531: PetscInt ngp;
532: PetscScalar gp_xi[GAUSS_POINTS][2];
533: PetscScalar gp_weight[GAUSS_POINTS];
534: PetscInt p,i,j,k;
535: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
536: PetscScalar J_p,tildeD[3];
537: PetscScalar B[3][U_DOFS*NODES_PER_EL];
540: /* define quadrature rule */
541: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
543: /* evaluate integral */
544: for (p = 0; p < ngp; p++) {
545: ConstructQ12D_GNi(gp_xi[p],GNi_p);
546: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
548: for (i = 0; i < NODES_PER_EL; i++) {
549: PetscScalar d_dx_i = GNx_p[0][i];
550: PetscScalar d_dy_i = GNx_p[1][i];
552: B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
553: B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
554: B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
555: }
558: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
559: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
560: tildeD[2] = gp_weight[p]*J_p*eta[p];
562: /* form Bt tildeD B */
563: /*
564: Ke_ij = Bt_ik . D_kl . B_lj
565: = B_ki . D_kl . B_lj
566: = B_ki . D_kk . B_kj
567: */
568: for (i = 0; i < 8; i++) {
569: for (j = 0; j < 8; j++) {
570: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
571: Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
572: }
573: }
574: }
575: }
576: }
578: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
579: {
580: PetscInt ngp;
581: PetscScalar gp_xi[GAUSS_POINTS][2];
582: PetscScalar gp_weight[GAUSS_POINTS];
583: PetscInt p,i,j,di;
584: PetscScalar Ni_p[NODES_PER_EL];
585: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
586: PetscScalar J_p,fac;
589: /* define quadrature rule */
590: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
592: /* evaluate integral */
593: for (p = 0; p < ngp; p++) {
594: ConstructQ12D_Ni(gp_xi[p],Ni_p);
595: ConstructQ12D_GNi(gp_xi[p],GNi_p);
596: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
597: fac = gp_weight[p]*J_p;
599: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
600: for (di = 0; di < NSD; di++) { /* u dofs */
601: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
602: PetscInt IJ;
603: /* Ke[4*u_idx+j] = Ke[4*u_idx+j] - GNx_p[di][i] * Ni_p[j] * fac; */
604: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);
606: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
607: }
608: }
609: }
610: }
611: }
613: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
614: {
615: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
616: PetscInt i,j;
617: PetscInt nr_g,nc_g;
619: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
620: FormGradientOperatorQ1(Ge,coords);
622: nr_g = U_DOFS*NODES_PER_EL;
623: nc_g = P_DOFS*NODES_PER_EL;
625: for (i = 0; i < nr_g; i++) {
626: for (j = 0; j < nc_g; j++) {
627: De[nr_g*j+i] = Ge[nc_g*i+j];
628: }
629: }
630: }
632: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
633: {
634: PetscInt ngp;
635: PetscScalar gp_xi[GAUSS_POINTS][2];
636: PetscScalar gp_weight[GAUSS_POINTS];
637: PetscInt p,i,j;
638: PetscScalar Ni_p[NODES_PER_EL];
639: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
640: PetscScalar J_p,fac,eta_avg;
643: /* define quadrature rule */
644: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
646: /* evaluate integral */
647: for (p = 0; p < ngp; p++) {
648: ConstructQ12D_Ni(gp_xi[p],Ni_p);
649: ConstructQ12D_GNi(gp_xi[p],GNi_p);
650: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
651: fac = gp_weight[p]*J_p;
653: for (i = 0; i < NODES_PER_EL; i++) {
654: for (j = 0; j < NODES_PER_EL; j++) {
655: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
656: }
657: }
658: }
660: /* scale */
661: eta_avg = 0.0;
662: for (p = 0; p < ngp; p++) {
663: eta_avg += eta[p];
664: }
665: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
666: fac = 1.0/eta_avg;
667: for (i = 0; i < NODES_PER_EL; i++) {
668: for (j = 0; j < NODES_PER_EL; j++) {
669: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
670: }
671: }
672: }
674: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
675: {
676: PetscInt ngp;
677: PetscScalar gp_xi[GAUSS_POINTS][2];
678: PetscScalar gp_weight[GAUSS_POINTS];
679: PetscInt p,i,j;
680: PetscScalar Ni_p[NODES_PER_EL];
681: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
682: PetscScalar J_p,fac,eta_avg;
685: /* define quadrature rule */
686: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
688: /* evaluate integral */
689: for (p = 0; p < ngp; p++) {
690: ConstructQ12D_Ni(gp_xi[p],Ni_p);
691: ConstructQ12D_GNi(gp_xi[p],GNi_p);
692: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
693: fac = gp_weight[p]*J_p;
695: for (i = 0; i < NODES_PER_EL; i++) {
696: for (j = 0; j < NODES_PER_EL; j++) {
697: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
698: }
699: }
700: }
702: /* scale */
703: eta_avg = 0.0;
704: for (p = 0; p < ngp; p++) {
705: eta_avg += eta[p];
706: }
707: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
708: fac = 1.0/eta_avg;
709: for (i = 0; i < NODES_PER_EL; i++) {
710: for (j = 0; j < NODES_PER_EL; j++) {
711: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
712: }
713: }
714: }
716: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
717: {
718: PetscInt ngp;
719: PetscScalar gp_xi[GAUSS_POINTS][2];
720: PetscScalar gp_weight[GAUSS_POINTS];
721: PetscInt p,i;
722: PetscScalar Ni_p[NODES_PER_EL];
723: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
724: PetscScalar J_p,fac;
727: /* define quadrature rule */
728: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
730: /* evaluate integral */
731: for (p = 0; p < ngp; p++) {
732: ConstructQ12D_Ni(gp_xi[p],Ni_p);
733: ConstructQ12D_GNi(gp_xi[p],GNi_p);
734: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
735: fac = gp_weight[p]*J_p;
737: for (i = 0; i < NODES_PER_EL; i++) {
738: Fe[NSD*i ] += fac*Ni_p[i]*fx[p];
739: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
740: }
741: }
742: }
746: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
747: {
749: /* get coords for the element */
750: el_coords[NSD*0+0] = _coords[ej][ei].x;el_coords[NSD*0+1] = _coords[ej][ei].y;
751: el_coords[NSD*1+0] = _coords[ej+1][ei].x;el_coords[NSD*1+1] = _coords[ej+1][ei].y;
752: el_coords[NSD*2+0] = _coords[ej+1][ei+1].x;el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
753: el_coords[NSD*3+0] = _coords[ej][ei+1].x;el_coords[NSD*3+1] = _coords[ej][ei+1].y;
754: return(0);
755: }
759: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
760: {
761: DM cda;
762: Vec coords;
763: DMDACoor2d **_coords;
764: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
765: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
766: PetscInt sex,sey,mx,my;
767: PetscInt ei,ej;
768: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
769: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
770: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
771: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
772: PetscScalar el_coords[NODES_PER_EL*NSD];
773: Vec local_properties;
774: GaussPointCoefficients **props;
775: PetscScalar *prop_eta;
776: PetscErrorCode ierr;
780: /* setup for coords */
781: DMDAGetCoordinateDA(stokes_da,&cda);
782: DMDAGetGhostedCoordinates(stokes_da,&coords);
783: DMDAVecGetArray(cda,coords,&_coords);
785: /* setup for coefficients */
786: DMCreateLocalVector(properties_da,&local_properties);
787: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
788: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
789: DMDAVecGetArray(properties_da,local_properties,&props);
791: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
792: for (ej = sey; ej < sey+my; ej++) {
793: for (ei = sex; ei < sex+mx; ei++) {
794: /* get coords for the element */
795: GetElementCoords(_coords,ei,ej,el_coords);
797: /* get coefficients for the element */
798: prop_eta = props[ej][ei].eta;
800: /* initialise element stiffness matrix */
801: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
802: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
803: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
804: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
806: /* form element stiffness matrix */
807: FormStressOperatorQ1(Ae,el_coords,prop_eta);
808: FormGradientOperatorQ1(Ge,el_coords);
809: FormDivergenceOperatorQ1(De,el_coords);
810: FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);
812: /* insert element matrix into global matrix */
813: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
814: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
815: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
816: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
817: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
818: }
819: }
820: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
821: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
823: DMDAVecRestoreArray(cda,coords,&_coords);
825: DMDAVecRestoreArray(properties_da,local_properties,&props);
826: VecDestroy(&local_properties);
827: return(0);
828: }
832: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
833: {
834: DM cda;
835: Vec coords;
836: DMDACoor2d **_coords;
837: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
838: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
839: PetscInt sex,sey,mx,my;
840: PetscInt ei,ej;
841: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
842: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
843: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
844: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
845: PetscScalar el_coords[NODES_PER_EL*NSD];
846: Vec local_properties;
847: GaussPointCoefficients **props;
848: PetscScalar *prop_eta;
849: PetscErrorCode ierr;
852: /* setup for coords */
853: DMDAGetCoordinateDA(stokes_da,&cda);
854: DMDAGetGhostedCoordinates(stokes_da,&coords);
855: DMDAVecGetArray(cda,coords,&_coords);
857: /* setup for coefficients */
858: DMCreateLocalVector(properties_da,&local_properties);
859: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
860: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
861: DMDAVecGetArray(properties_da,local_properties,&props);
863: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
864: for (ej = sey; ej < sey+my; ej++) {
865: for (ei = sex; ei < sex+mx; ei++) {
866: /* get coords for the element */
867: GetElementCoords(_coords,ei,ej,el_coords);
869: /* get coefficients for the element */
870: prop_eta = props[ej][ei].eta;
872: /* initialise element stiffness matrix */
873: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
874: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
875: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
876: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
879: /* form element stiffness matrix */
880: FormStressOperatorQ1(Ae,el_coords,prop_eta);
881: FormGradientOperatorQ1(Ge,el_coords);
882: /* FormDivergenceOperatorQ1( De, el_coords ); */
883: FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);
885: /* insert element matrix into global matrix */
886: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
887: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
888: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
889: /* MatSetValuesStencil( A, NODES_PER_EL*P_DOFS,p_eqn, NODES_PER_EL*U_DOFS,u_eqn, De, ADD_VALUES ); */
890: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
891: }
892: }
893: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
894: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
896: DMDAVecRestoreArray(cda,coords,&_coords);
898: DMDAVecRestoreArray(properties_da,local_properties,&props);
899: VecDestroy(&local_properties);
900: return(0);
901: }
905: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
906: {
907: PetscInt n;
910: for (n = 0; n < 4; n++) {
911: fields_F[ u_eqn[2*n ].j ][ u_eqn[2*n ].i ].u_dof = fields_F[ u_eqn[2*n ].j ][ u_eqn[2*n ].i ].u_dof+Fe_u[2*n ];
912: fields_F[ u_eqn[2*n+1].j ][ u_eqn[2*n+1].i ].v_dof = fields_F[ u_eqn[2*n+1].j ][ u_eqn[2*n+1].i ].v_dof+Fe_u[2*n+1];
913: fields_F[ p_eqn[n].j ][ p_eqn[n].i ].p_dof = fields_F[ p_eqn[n].j ][ p_eqn[n].i ].p_dof+Fe_p[n ];
914: }
915: return(0);
916: }
920: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
921: {
922: DM cda;
923: Vec coords;
924: DMDACoor2d **_coords;
925: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
926: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
927: PetscInt sex,sey,mx,my;
928: PetscInt ei,ej;
929: PetscScalar Fe[NODES_PER_EL*U_DOFS];
930: PetscScalar He[NODES_PER_EL*P_DOFS];
931: PetscScalar el_coords[NODES_PER_EL*NSD];
932: Vec local_properties;
933: GaussPointCoefficients **props;
934: PetscScalar *prop_fx,*prop_fy;
935: Vec local_F;
936: StokesDOF **ff;
937: PetscErrorCode ierr;
940: /* setup for coords */
941: DMDAGetCoordinateDA(stokes_da,&cda);
942: DMDAGetGhostedCoordinates(stokes_da,&coords);
943: DMDAVecGetArray(cda,coords,&_coords);
945: /* setup for coefficients */
946: DMGetLocalVector(properties_da,&local_properties);
947: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
948: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
949: DMDAVecGetArray(properties_da,local_properties,&props);
951: /* get acces to the vector */
952: DMGetLocalVector(stokes_da,&local_F);
953: VecZeroEntries(local_F);
954: DMDAVecGetArray(stokes_da,local_F,&ff);
957: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
958: for (ej = sey; ej < sey+my; ej++) {
959: for (ei = sex; ei < sex+mx; ei++) {
960: /* get coords for the element */
961: GetElementCoords(_coords,ei,ej,el_coords);
963: /* get coefficients for the element */
964: prop_fx = props[ej][ei].fx;
965: prop_fy = props[ej][ei].fy;
967: /* initialise element stiffness matrix */
968: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
969: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
972: /* form element stiffness matrix */
973: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
975: /* insert element matrix into global matrix */
976: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
978: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
979: }
980: }
982: DMDAVecRestoreArray(stokes_da,local_F,&ff);
983: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
984: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
985: DMRestoreLocalVector(stokes_da,&local_F);
988: DMDAVecRestoreArray(cda,coords,&_coords);
990: DMDAVecRestoreArray(properties_da,local_properties,&props);
991: DMRestoreLocalVector(properties_da,&local_properties);
992: return(0);
993: }
997: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,
998: PetscInt mx,PetscInt my,
999: DM *_da,Vec *_X)
1000: {
1001: DM da,cda;
1002: Vec X,local_X;
1003: StokesDOF **_stokes;
1004: Vec coords;
1005: DMDACoor2d **_coords;
1006: PetscInt si,sj,ei,ej,i,j;
1010: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1011: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,3,1,PETSC_NULL,PETSC_NULL,&da);
1012: DMDASetFieldName(da,0,"anlytic_Vx");
1013: DMDASetFieldName(da,1,"anlytic_Vy");
1014: DMDASetFieldName(da,2,"analytic_P");
1017: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,PETSC_NULL,PETSC_NULL);
1020: DMDAGetGhostedCoordinates(da,&coords);
1021: DMDAGetCoordinateDA(da,&cda);
1022: DMDAVecGetArray(cda,coords,&_coords);
1024: DMCreateGlobalVector(da,&X);
1025: DMCreateLocalVector(da,&local_X);
1026: DMDAVecGetArray(da,local_X,&_stokes);
1028: DMDAGetGhostCorners(da,&si,&sj,0,&ei,&ej,0);
1029: for (j = sj; j < sj+ej; j++) {
1030: for (i = si; i < si+ei; i++) {
1031: double pos[2],pressure,vel[2],total_stress[3],strain_rate[3];
1033: pos[0] = PetscRealPart(_coords[j][i].x);
1034: pos[1] = PetscRealPart(_coords[j][i].y);
1036: evaluate_solCx(pos,eta0,eta1,xc,nz,vel,&pressure,total_stress,strain_rate);
1038: _stokes[j][i].u_dof = vel[0];
1039: _stokes[j][i].v_dof = vel[1];
1040: _stokes[j][i].p_dof = pressure;
1041: }
1042: }
1043: DMDAVecRestoreArray(da,local_X,&_stokes);
1044: DMDAVecRestoreArray(cda,coords,&_coords);
1046: DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1047: DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);
1049: VecDestroy(&local_X);
1051: *_da = da;
1052: *_X = X;
1054: return(0);
1055: }
1059: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
1060: {
1062: /* get the nodal fields */
1063: nodal_fields[0].u_dof = fields[ej ][ei ].u_dof;nodal_fields[0].v_dof = fields[ej ][ei ].v_dof;nodal_fields[0].p_dof = fields[ej ][ei ].p_dof;
1064: nodal_fields[1].u_dof = fields[ej+1][ei ].u_dof;nodal_fields[1].v_dof = fields[ej+1][ei ].v_dof;nodal_fields[1].p_dof = fields[ej+1][ei ].p_dof;
1065: nodal_fields[2].u_dof = fields[ej+1][ei+1].u_dof;nodal_fields[2].v_dof = fields[ej+1][ei+1].v_dof;nodal_fields[2].p_dof = fields[ej+1][ei+1].p_dof;
1066: nodal_fields[3].u_dof = fields[ej ][ei+1].u_dof;nodal_fields[3].v_dof = fields[ej ][ei+1].v_dof;nodal_fields[3].p_dof = fields[ej ][ei+1].p_dof;
1067: return(0);
1068: }
1072: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
1073: {
1074: DM cda;
1075: Vec coords,X_analytic_local,X_local;
1076: DMDACoor2d **_coords;
1077: PetscInt sex,sey,mx,my;
1078: PetscInt ei,ej;
1079: PetscScalar el_coords[NODES_PER_EL*NSD];
1080: StokesDOF **stokes_analytic,**stokes;
1081: StokesDOF stokes_analytic_e[4],stokes_e[4];
1083: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1084: PetscScalar Ni_p[NODES_PER_EL];
1085: PetscInt ngp;
1086: PetscScalar gp_xi[GAUSS_POINTS][2];
1087: PetscScalar gp_weight[GAUSS_POINTS];
1088: PetscInt p,i;
1089: PetscScalar J_p,fac;
1090: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1091: PetscInt M;
1092: PetscReal xymin[2],xymax[2];
1096: /* define quadrature rule */
1097: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1099: /* setup for coords */
1100: DMDAGetCoordinateDA(stokes_da,&cda);
1101: DMDAGetGhostedCoordinates(stokes_da,&coords);
1102: DMDAVecGetArray(cda,coords,&_coords);
1104: /* setup for analytic */
1105: DMCreateLocalVector(stokes_da,&X_analytic_local);
1106: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1107: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1108: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1110: /* setup for solution */
1111: DMCreateLocalVector(stokes_da,&X_local);
1112: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1113: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1114: DMDAVecGetArray(stokes_da,X_local,&stokes);
1116: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1117: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1119: h = (xymax[0]-xymin[0])/((double)M);
1121: tp_L2 = tu_L2 = tu_H1 = 0.0;
1123: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
1124: for (ej = sey; ej < sey+my; ej++) {
1125: for (ei = sex; ei < sex+mx; ei++) {
1126: /* get coords for the element */
1127: GetElementCoords(_coords,ei,ej,el_coords);
1128: StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1129: StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);
1131: /* evaluate integral */
1132: p_e_L2 = 0.0;
1133: u_e_L2 = 0.0;
1134: u_e_H1 = 0.0;
1135: for (p = 0; p < ngp; p++) {
1136: ConstructQ12D_Ni(gp_xi[p],Ni_p);
1137: ConstructQ12D_GNi(gp_xi[p],GNi_p);
1138: ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1139: fac = gp_weight[p]*J_p;
1141: for (i = 0; i < NODES_PER_EL; i++) {
1142: PetscScalar u_error,v_error;
1144: p_e_L2 = p_e_L2+fac*Ni_p[i]*(stokes_e[i].p_dof-stokes_analytic_e[i].p_dof)*(stokes_e[i].p_dof-stokes_analytic_e[i].p_dof);
1146: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1147: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1148: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);
1150: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1151: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error /* du/dy */
1152: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1153: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error); /* dv/dy */
1154: }
1155: }
1157: tp_L2 += p_e_L2;
1158: tu_L2 += u_e_L2;
1159: tu_H1 += u_e_H1;
1160: }
1161: }
1162: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1163: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1164: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1165: p_L2 = PetscSqrtScalar(p_L2);
1166: u_L2 = PetscSqrtScalar(u_L2);
1167: u_H1 = PetscSqrtScalar(u_H1);
1169: PetscPrintf(PETSC_COMM_WORLD,"%1.4e %1.4e %1.4e %1.4e \n",PetscRealPart(h),PetscRealPart(p_L2),PetscRealPart(u_L2),PetscRealPart(u_H1));
1172: DMDAVecRestoreArray(cda,coords,&_coords);
1174: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1175: VecDestroy(&X_analytic_local);
1176: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1177: VecDestroy(&X_local);
1178: return(0);
1179: }
1183: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1184: {
1185: DM da_Stokes,da_prop;
1186: PetscInt u_dof,p_dof,dof,stencil_width;
1187: Mat A,B;
1188: PetscInt mxl,myl;
1189: DM prop_cda,vel_cda;
1190: Vec prop_coords,vel_coords;
1191: PetscInt si,sj,nx,ny,i,j,p;
1192: Vec f,X;
1193: PetscInt prop_dof,prop_stencil_width;
1194: Vec properties,l_properties;
1195: PetscReal dx,dy;
1196: PetscInt M,N;
1197: DMDACoor2d **_prop_coords,**_vel_coords;
1198: GaussPointCoefficients **element_props;
1199: PetscInt its;
1200: KSP ksp_S;
1201: PetscInt coefficient_structure = 0;
1202: PetscInt cpu_x,cpu_y,*lx = PETSC_NULL,*ly = PETSC_NULL;
1203: PetscBool use_gp_coords = PETSC_FALSE,set;
1204: char filename[PETSC_MAX_PATH_LEN];
1205: PetscErrorCode ierr;
1208: /* Generate the da for velocity and pressure */
1209: /*
1210: We use Q1 elements for the temperature.
1211: FEM has a 9-point stencil (BOX) or connectivity pattern
1212: Num nodes in each direction is mx+1, my+1
1213: */
1214: u_dof = U_DOFS; /* Vx, Vy - velocities */
1215: p_dof = P_DOFS; /* p - pressure */
1216: dof = u_dof+p_dof;
1217: stencil_width = 1;
1218: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1219: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,PETSC_NULL,PETSC_NULL,&da_Stokes);
1220: DMDASetFieldName(da_Stokes,0,"Vx");
1221: DMDASetFieldName(da_Stokes,1,"Vy");
1222: DMDASetFieldName(da_Stokes,2,"P");
1224: /* unit box [0,1] x [0,1] */
1225: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,PETSC_NULL,PETSC_NULL);
1228: /* Generate element properties, we will assume all material properties are constant over the element */
1229: /* local number of elements */
1230: DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,PETSC_NULL);
1232: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! // */
1233: DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1234: DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);
1236: prop_dof = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1237: prop_stencil_width = 0;
1238: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1239: mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
1240: PetscFree(lx);
1241: PetscFree(ly);
1243: /* define centroid positions */
1244: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1245: dx = 1.0/((PetscReal)(M));
1246: dy = 1.0/((PetscReal)(N));
1248: DMDASetUniformCoordinates(da_prop,0.0+0.5*dx,1.0-0.5*dx,0.0+0.5*dy,1.0-0.5*dy,PETSC_NULL,PETSC_NULL);
1250: /* define coefficients */
1251: PetscOptionsGetInt(PETSC_NULL,"-c_str",&coefficient_structure,PETSC_NULL);
1252: /* PetscPrintf( PETSC_COMM_WORLD, "Using coeficient structure %D \n", coefficient_structure ); */
1254: DMCreateGlobalVector(da_prop,&properties);
1255: DMCreateLocalVector(da_prop,&l_properties);
1256: DMDAVecGetArray(da_prop,l_properties,&element_props);
1258: DMDAGetCoordinateDA(da_prop,&prop_cda);
1259: DMDAGetGhostedCoordinates(da_prop,&prop_coords);
1260: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
1262: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
1264: DMDAGetCoordinateDA(da_Stokes,&vel_cda);
1265: DMDAGetGhostedCoordinates(da_Stokes,&vel_coords);
1266: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1269: /* interpolate the coordinates */
1270: for (j = sj; j < sj+ny; j++) {
1271: for (i = si; i < si+nx; i++) {
1272: PetscInt ngp;
1273: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1274: PetscScalar el_coords[8];
1276: GetElementCoords(_vel_coords,i,j,el_coords);
1277: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1279: for (p = 0; p < GAUSS_POINTS; p++) {
1280: PetscScalar gp_x,gp_y;
1281: PetscInt n;
1282: PetscScalar xi_p[2],Ni_p[4];
1284: xi_p[0] = gp_xi[p][0];
1285: xi_p[1] = gp_xi[p][1];
1286: ConstructQ12D_Ni(xi_p,Ni_p);
1288: gp_x = 0.0;
1289: gp_y = 0.0;
1290: for (n = 0; n < NODES_PER_EL; n++) {
1291: gp_x = gp_x+Ni_p[n]*el_coords[2*n ];
1292: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1293: }
1294: element_props[j][i].gp_coords[2*p ] = gp_x;
1295: element_props[j][i].gp_coords[2*p+1] = gp_y;
1296: }
1297: }
1298: }
1300: /* define the coefficients */
1301: PetscOptionsGetBool(PETSC_NULL,"-use_gp_coords",&use_gp_coords,0);
1303: for (j = sj; j < sj+ny; j++) {
1304: for (i = si; i < si+nx; i++) {
1305: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1306: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1307: PetscReal coord_x,coord_y;
1309: if (coefficient_structure == 0) {
1310: PetscReal opts_eta0,opts_eta1,opts_xc;
1311: PetscInt opts_nz;
1313: opts_eta0 = 1.0;
1314: opts_eta1 = 1.0;
1315: opts_xc = 0.5;
1316: opts_nz = 1;
1318: PetscOptionsGetReal(PETSC_NULL,"-solcx_eta0",&opts_eta0,0);
1319: PetscOptionsGetReal(PETSC_NULL,"-solcx_eta1",&opts_eta1,0);
1320: PetscOptionsGetReal(PETSC_NULL,"-solcx_xc",&opts_xc,0);
1321: PetscOptionsGetInt(PETSC_NULL,"-solcx_nz",&opts_nz,0);
1323: for (p = 0; p < GAUSS_POINTS; p++) {
1324: coord_x = centroid_x;
1325: coord_y = centroid_y;
1326: if (use_gp_coords == PETSC_TRUE) {
1327: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1328: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1329: }
1332: element_props[j][i].eta[p] = opts_eta0;
1333: if (coord_x > opts_xc) {
1334: element_props[j][i].eta[p] = opts_eta1;
1335: }
1337: element_props[j][i].fx[p] = 0.0;
1338: element_props[j][i].fy[p] = sin((double)opts_nz*PETSC_PI*coord_y)*cos(1.0*PETSC_PI*coord_x);
1339: }
1340: } else if (coefficient_structure == 1) { /* square sinker */
1341: PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;
1343: opts_eta0 = 1.0;
1344: opts_eta1 = 1.0;
1345: opts_dx = 0.50;
1346: opts_dy = 0.50;
1348: PetscOptionsGetReal(PETSC_NULL,"-sinker_eta0",&opts_eta0,0);
1349: PetscOptionsGetReal(PETSC_NULL,"-sinker_eta1",&opts_eta1,0);
1350: PetscOptionsGetReal(PETSC_NULL,"-sinker_dx",&opts_dx,0);
1351: PetscOptionsGetReal(PETSC_NULL,"-sinker_dy",&opts_dy,0);
1354: for (p = 0; p < GAUSS_POINTS; p++) {
1355: coord_x = centroid_x;
1356: coord_y = centroid_y;
1357: if (use_gp_coords == PETSC_TRUE) {
1358: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1359: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1360: }
1362: element_props[j][i].eta[p] = opts_eta0;
1363: element_props[j][i].fx[p] = 0.0;
1364: element_props[j][i].fy[p] = 0.0;
1366: if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1367: if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1368: element_props[j][i].eta[p] = opts_eta1;
1369: element_props[j][i].fx[p] = 0.0;
1370: element_props[j][i].fy[p] = -1.0;
1371: }
1372: }
1373: }
1374: } else if (coefficient_structure == 2) { /* circular sinker */
1375: PetscReal opts_eta0,opts_eta1,opts_r,radius2;
1377: opts_eta0 = 1.0;
1378: opts_eta1 = 1.0;
1379: opts_r = 0.25;
1381: PetscOptionsGetReal(PETSC_NULL,"-sinker_eta0",&opts_eta0,0);
1382: PetscOptionsGetReal(PETSC_NULL,"-sinker_eta1",&opts_eta1,0);
1383: PetscOptionsGetReal(PETSC_NULL,"-sinker_r",&opts_r,0);
1385: for (p = 0; p < GAUSS_POINTS; p++) {
1386: coord_x = centroid_x;
1387: coord_y = centroid_y;
1388: if (use_gp_coords == PETSC_TRUE) {
1389: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1390: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1391: }
1393: element_props[j][i].eta[p] = opts_eta0;
1394: element_props[j][i].fx[p] = 0.0;
1395: element_props[j][i].fy[p] = 0.0;
1397: radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1398: if (radius2 < opts_r*opts_r) {
1399: element_props[j][i].eta[p] = opts_eta1;
1400: element_props[j][i].fx[p] = 0.0;
1401: element_props[j][i].fy[p] = -1.0;
1402: }
1403: }
1404: } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1405: PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;
1407: opts_eta0 = 1.0;
1408: opts_eta1 = 1.0;
1409: opts_r = 0.25;
1410: opts_c0x = 0.35; /* circle center */
1411: opts_c0y = 0.35;
1412: opts_s0x = 0.7; /* square center */
1413: opts_s0y = 0.7;
1414: opts_dx = 0.25;
1415: opts_dy = 0.25;
1416: opts_phi = 25;
1418: PetscOptionsGetReal(PETSC_NULL,"-sinker_eta0",&opts_eta0,0);
1419: PetscOptionsGetReal(PETSC_NULL,"-sinker_eta1",&opts_eta1,0);
1420: PetscOptionsGetReal(PETSC_NULL,"-sinker_r",&opts_r,0);
1421: PetscOptionsGetReal(PETSC_NULL,"-sinker_c0x",&opts_c0x,0);
1422: PetscOptionsGetReal(PETSC_NULL,"-sinker_c0y",&opts_c0y,0);
1423: PetscOptionsGetReal(PETSC_NULL,"-sinker_s0x",&opts_s0x,0);
1424: PetscOptionsGetReal(PETSC_NULL,"-sinker_s0y",&opts_s0y,0);
1425: PetscOptionsGetReal(PETSC_NULL,"-sinker_dx",&opts_dx,0);
1426: PetscOptionsGetReal(PETSC_NULL,"-sinker_dy",&opts_dy,0);
1427: PetscOptionsGetReal(PETSC_NULL,"-sinker_phi",&opts_phi,0);
1428: opts_phi *= PETSC_PI / 180;
1430: for (p = 0; p < GAUSS_POINTS; p++) {
1431: coord_x = centroid_x;
1432: coord_y = centroid_y;
1433: if (use_gp_coords == PETSC_TRUE) {
1434: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1435: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1436: }
1438: element_props[j][i].eta[p] = opts_eta0;
1439: element_props[j][i].fx[p] = 0.0;
1440: element_props[j][i].fy[p] = -0.2;
1442: radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1443: if (radius2 < opts_r*opts_r
1444: || ( PetscAbs(+(coord_x - opts_s0x)*cos(opts_phi) + (coord_y - opts_s0y)*sin(opts_phi)) < opts_dx/2
1445: && PetscAbs(-(coord_x - opts_s0x)*sin(opts_phi) + (coord_y - opts_s0y)*cos(opts_phi)) < opts_dy/2)) {
1446: element_props[j][i].eta[p] = opts_eta1;
1447: element_props[j][i].fx[p] = 0.0;
1448: element_props[j][i].fy[p] = -1.0;
1449: }
1450: }
1451: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1452: }
1453: }
1454: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
1456: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1458: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1459: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1460: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
1463: DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1464: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for Stokes eqn.","properties");
1467: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1468: DMCreateMatrix(da_Stokes,MATAIJ,&A);
1469: DMCreateMatrix(da_Stokes,MATAIJ,&B);
1470: DMCreateGlobalVector(da_Stokes,&f);
1471: DMCreateGlobalVector(da_Stokes,&X);
1473: /* assemble A11 */
1474: MatZeroEntries(A);
1475: MatZeroEntries(B);
1476: VecZeroEntries(f);
1478: AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1479: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1480: /* build force vector */
1481: AssembleF_Stokes(f,da_Stokes,da_prop,properties);
1483: DMDABCApplyFreeSlip(da_Stokes,A,f);
1484: DMDABCApplyFreeSlip(da_Stokes,B,PETSC_NULL);
1486: /* SOLVE */
1487: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1488: KSPSetOptionsPrefix(ksp_S,"stokes_");
1489: KSPSetOperators(ksp_S,A,B,SAME_NONZERO_PATTERN);
1490: KSPSetDM(ksp_S,da_Stokes);
1491: KSPSetDMActive(ksp_S,PETSC_FALSE);
1492: KSPSetFromOptions(ksp_S);
1493: {
1494: PC pc;
1495: const PetscInt ufields[] = {0,1},pfields[1] = {2};
1496: KSPGetPC(ksp_S,&pc);
1497: PCFieldSplitSetBlockSize(pc,3);
1498: PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1499: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1500: }
1502: KSPSolve(ksp_S,f,X);
1504: PetscOptionsGetString(PETSC_NULL,"-o",filename,sizeof filename,&set);
1505: if (set) {
1506: char *ext;
1507: PetscViewer viewer;
1508: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1509: PetscStrrchr(filename,'.',&ext);
1510: if (!strcmp("vts",ext)) {
1511: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1512: } else {
1513: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1514: }
1515: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1516: PetscViewerFileSetName(viewer,filename);
1517: VecView(X,viewer);
1518: PetscViewerDestroy(&viewer);
1519: }
1520: DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
1522: KSPGetIterationNumber(ksp_S,&its);
1524: /*
1525: {
1526: Vec x;
1527: PetscInt L;
1528: VecDuplicate(f,&x);
1529: MatMult(A,X, x );
1530: VecAXPY( x, -1.0, f );
1531: VecNorm( x, NORM_2, &nrm );
1532: PetscPrintf( PETSC_COMM_WORLD, "Its. %1.4d, norm |AX-f| = %1.5e \n", its, nrm );
1533: VecDestroy(&x);
1535: VecNorm( X, NORM_2, &nrm );
1536: VecGetSize( X, &L );
1537: PetscPrintf( PETSC_COMM_WORLD, " norm |X|/sqrt{N} = %1.5e \n", nrm/sqrt( (PetscScalar)L ) );
1538: }
1539: */
1541: if (coefficient_structure == 0) {
1542: PetscReal opts_eta0,opts_eta1,opts_xc;
1543: PetscInt opts_nz,N;
1544: DM da_Stokes_analytic;
1545: Vec X_analytic;
1546: PetscReal nrm1[3],nrm2[3],nrmI[3];
1548: opts_eta0 = 1.0;
1549: opts_eta1 = 1.0;
1550: opts_xc = 0.5;
1551: opts_nz = 1;
1553: PetscOptionsGetReal(PETSC_NULL,"-solcx_eta0",&opts_eta0,0);
1554: PetscOptionsGetReal(PETSC_NULL,"-solcx_eta1",&opts_eta1,0);
1555: PetscOptionsGetReal(PETSC_NULL,"-solcx_xc",&opts_xc,0);
1556: PetscOptionsGetInt(PETSC_NULL,"-solcx_nz",&opts_nz,0);
1559: DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1560: DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
1562: DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);
1565: VecAXPY(X_analytic,-1.0,X);
1566: VecGetSize(X_analytic,&N);
1567: N = N/3;
1569: VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1570: VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1571: VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);
1573: VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1574: VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1575: VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);
1577: VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1578: VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1579: VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);
1581: /*
1582: PetscPrintf( PETSC_COMM_WORLD, "%1.4e %1.4e %1.4e %1.4e %1.4e %1.4e %1.4e %1.4e %1.4e %1.4e\n",
1583: 1.0/mx,
1584: nrm1[0]/(double)N, sqrt(nrm2[0]/(double)N), nrmI[0],
1585: nrm1[1]/(double)N, sqrt(nrm2[1]/(double)N), nrmI[1],
1586: nrm1[2]/(double)N, sqrt(nrm2[2]/(double)N), nrmI[2] );
1587: */
1588: DMDestroy(&da_Stokes_analytic);
1589: VecDestroy(&X_analytic);
1590: }
1593: KSPDestroy(&ksp_S);
1594: VecDestroy(&X);
1595: VecDestroy(&f);
1596: MatDestroy(&A);
1597: MatDestroy(&B);
1599: DMDestroy(&da_Stokes);
1600: DMDestroy(&da_prop);
1602: VecDestroy(&properties);
1603: VecDestroy(&l_properties);
1605: return(0);
1606: }
1610: int main(int argc,char **args)
1611: {
1613: PetscInt mx,my;
1615: PetscInitialize(&argc,&args,(char *)0,help);
1617: mx = my = 10;
1618: PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,PETSC_NULL);
1619: PetscOptionsGetInt(PETSC_NULL,"-my",&my,PETSC_NULL);
1621: solve_stokes_2d_coupled(mx,my);
1623: PetscFinalize();
1624: return 0;
1625: }
1627: /* -------------------------- helpers for boundary conditions -------------------------------- */
1631: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1632: {
1633: DM cda;
1634: Vec coords;
1635: PetscInt si,sj,nx,ny,i,j;
1636: PetscInt M,N;
1637: DMDACoor2d **_coords;
1638: PetscInt *g_idx;
1639: PetscInt *bc_global_ids;
1640: PetscScalar *bc_vals;
1641: PetscInt nbcs;
1642: PetscInt n_dofs;
1646: /* enforce bc's */
1647: DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);
1649: DMDAGetCoordinateDA(da,&cda);
1650: DMDAGetGhostedCoordinates(da,&coords);
1651: DMDAVecGetArray(cda,coords,&_coords);
1652: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1653: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1655: /* /// */
1657: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1658: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1660: /* init the entries to -1 so VecSetValues will ignore them */
1661: for (i = 0; i < ny*n_dofs; i++) {
1662: bc_global_ids[i] = -1;
1663: }
1665: i = nx-1;
1666: for (j = 0; j < ny; j++) {
1667: PetscInt local_id;
1669: local_id = i+j*nx;
1671: bc_global_ids[j] = g_idx[ n_dofs*local_id+d_idx ];
1673: bc_vals[j] = bc_val;
1674: }
1675: nbcs = 0;
1676: if ((si+nx) == (M)) {
1677: nbcs = ny;
1678: }
1680: if (b != PETSC_NULL) {
1681: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1682: VecAssemblyBegin(b);
1683: VecAssemblyEnd(b);
1684: }
1685: if (A != PETSC_NULL) {
1686: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1687: }
1690: PetscFree(bc_vals);
1691: PetscFree(bc_global_ids);
1693: DMDAVecRestoreArray(cda,coords,&_coords);
1694: return(0);
1695: }
1699: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1700: {
1701: DM cda;
1702: Vec coords;
1703: PetscInt si,sj,nx,ny,i,j;
1704: PetscInt M,N;
1705: DMDACoor2d **_coords;
1706: PetscInt *g_idx;
1707: PetscInt *bc_global_ids;
1708: PetscScalar *bc_vals;
1709: PetscInt nbcs;
1710: PetscInt n_dofs;
1714: /* enforce bc's */
1715: DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);
1717: DMDAGetCoordinateDA(da,&cda);
1718: DMDAGetGhostedCoordinates(da,&coords);
1719: DMDAVecGetArray(cda,coords,&_coords);
1720: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1721: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1723: /* /// */
1725: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1726: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1728: /* init the entries to -1 so VecSetValues will ignore them */
1729: for (i = 0; i < ny*n_dofs; i++) {
1730: bc_global_ids[i] = -1;
1731: }
1733: i = 0;
1734: for (j = 0; j < ny; j++) {
1735: PetscInt local_id;
1737: local_id = i+j*nx;
1739: bc_global_ids[j] = g_idx[ n_dofs*local_id+d_idx ];
1741: bc_vals[j] = bc_val;
1742: }
1743: nbcs = 0;
1744: if (si == 0) {
1745: nbcs = ny;
1746: }
1748: if (b != PETSC_NULL) {
1749: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1750: VecAssemblyBegin(b);
1751: VecAssemblyEnd(b);
1752: }
1753: if (A != PETSC_NULL) {
1754: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1755: }
1758: PetscFree(bc_vals);
1759: PetscFree(bc_global_ids);
1761: DMDAVecRestoreArray(cda,coords,&_coords);
1762: return(0);
1763: }
1767: static PetscErrorCode BCApply_NORTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1768: {
1769: DM cda;
1770: Vec coords;
1771: PetscInt si,sj,nx,ny,i,j;
1772: PetscInt M,N;
1773: DMDACoor2d **_coords;
1774: PetscInt *g_idx;
1775: PetscInt *bc_global_ids;
1776: PetscScalar *bc_vals;
1777: PetscInt nbcs;
1778: PetscInt n_dofs;
1782: /* enforce bc's */
1783: DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);
1785: DMDAGetCoordinateDA(da,&cda);
1786: DMDAGetGhostedCoordinates(da,&coords);
1787: DMDAVecGetArray(cda,coords,&_coords);
1788: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1789: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1791: /* /// */
1793: PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1794: PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);
1796: /* init the entries to -1 so VecSetValues will ignore them */
1797: for (i = 0; i < nx; i++) {
1798: bc_global_ids[i] = -1;
1799: }
1801: j = ny-1;
1802: for (i = 0; i < nx; i++) {
1803: PetscInt local_id;
1805: local_id = i+j*nx;
1807: bc_global_ids[i] = g_idx[ n_dofs*local_id+d_idx ];
1809: bc_vals[i] = bc_val;
1810: }
1811: nbcs = 0;
1812: if ((sj+ny) == (N)) {
1813: nbcs = nx;
1814: }
1816: if (b != PETSC_NULL) {
1817: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1818: VecAssemblyBegin(b);
1819: VecAssemblyEnd(b);
1820: }
1821: if (A != PETSC_NULL) {
1822: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1823: }
1826: PetscFree(bc_vals);
1827: PetscFree(bc_global_ids);
1829: DMDAVecRestoreArray(cda,coords,&_coords);
1830: return(0);
1831: }
1835: static PetscErrorCode BCApply_SOUTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1836: {
1837: DM cda;
1838: Vec coords;
1839: PetscInt si,sj,nx,ny,i,j;
1840: PetscInt M,N;
1841: DMDACoor2d **_coords;
1842: PetscInt *g_idx;
1843: PetscInt *bc_global_ids;
1844: PetscScalar *bc_vals;
1845: PetscInt nbcs;
1846: PetscInt n_dofs;
1850: /* enforce bc's */
1851: DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);
1853: DMDAGetCoordinateDA(da,&cda);
1854: DMDAGetGhostedCoordinates(da,&coords);
1855: DMDAVecGetArray(cda,coords,&_coords);
1856: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1857: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1859: /* /// */
1861: PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1862: PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);
1864: /* init the entries to -1 so VecSetValues will ignore them */
1865: for (i = 0; i < nx; i++) {
1866: bc_global_ids[i] = -1;
1867: }
1869: j = 0;
1870: for (i = 0; i < nx; i++) {
1871: PetscInt local_id;
1873: local_id = i+j*nx;
1875: bc_global_ids[i] = g_idx[ n_dofs*local_id+d_idx ];
1877: bc_vals[i] = bc_val;
1878: }
1879: nbcs = 0;
1880: if (sj == 0) {
1881: nbcs = nx;
1882: }
1885: if (b != PETSC_NULL) {
1886: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1887: VecAssemblyBegin(b);
1888: VecAssemblyEnd(b);
1889: }
1890: if (A != PETSC_NULL) {
1891: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1892: }
1895: PetscFree(bc_vals);
1896: PetscFree(bc_global_ids);
1898: DMDAVecRestoreArray(cda,coords,&_coords);
1899: return(0);
1900: }
1904: /*
1905: Free slip sides.
1906: */
1909: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1910: {
1914: BCApply_NORTH(da_Stokes,1,0.0,A,f);
1915: BCApply_EAST(da_Stokes,0,0.0,A,f);
1916: BCApply_SOUTH(da_Stokes,1,0.0,A,f);
1917: BCApply_WEST(da_Stokes,0,0.0,A,f);
1918: return(0);
1919: }