Actual source code: ex43.c
petsc-3.4.5 2014-06-29
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 <petscdmda.h>
36: /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
37: #include "ex43-solcx.h"
39: static PetscErrorCode DMDABCApplyFreeSlip(DM,Mat,Vec);
42: #define NSD 2 /* number of spatial dimensions */
43: #define NODES_PER_EL 4 /* nodes per element */
44: #define U_DOFS 2 /* degrees of freedom per velocity node */
45: #define P_DOFS 1 /* degrees of freedom per pressure node */
46: #define GAUSS_POINTS 4
48: /* cell based evaluation */
49: typedef struct {
50: PetscScalar eta,fx,fy;
51: } Coefficients;
53: /* Gauss point based evaluation 8+4+4+4 = 20 */
54: typedef struct {
55: PetscScalar gp_coords[2*GAUSS_POINTS];
56: PetscScalar eta[GAUSS_POINTS];
57: PetscScalar fx[GAUSS_POINTS];
58: PetscScalar fy[GAUSS_POINTS];
59: } GaussPointCoefficients;
61: typedef struct {
62: PetscScalar u_dof;
63: PetscScalar v_dof;
64: PetscScalar p_dof;
65: } StokesDOF;
68: /*
70: D = [ 2.eta 0 0 ]
71: [ 0 2.eta 0 ]
72: [ 0 0 eta ]
74: B = [ d_dx 0 ]
75: [ 0 d_dy ]
76: [ d_dy d_dx ]
78: */
80: /* FEM routines */
81: /*
82: Element: Local basis function ordering
83: 1-----2
84: | |
85: | |
86: 0-----3
87: */
88: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
89: {
90: PetscScalar xi = _xi[0];
91: PetscScalar eta = _xi[1];
93: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
94: Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
95: Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
96: Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
97: }
99: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
100: {
101: PetscScalar xi = _xi[0];
102: PetscScalar eta = _xi[1];
104: GNi[0][0] = -0.25*(1.0-eta);
105: GNi[0][1] = -0.25*(1.0+eta);
106: GNi[0][2] = 0.25*(1.0+eta);
107: GNi[0][3] = 0.25*(1.0-eta);
109: GNi[1][0] = -0.25*(1.0-xi);
110: GNi[1][1] = 0.25*(1.0-xi);
111: GNi[1][2] = 0.25*(1.0+xi);
112: GNi[1][3] = -0.25*(1.0+xi);
113: }
115: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
116: {
117: PetscScalar J00,J01,J10,J11,J;
118: PetscScalar iJ00,iJ01,iJ10,iJ11;
119: PetscInt i;
121: J00 = J01 = J10 = J11 = 0.0;
122: for (i = 0; i < NODES_PER_EL; i++) {
123: PetscScalar cx = coords[2*i+0];
124: PetscScalar cy = coords[2*i+1];
126: J00 = J00+GNi[0][i]*cx; /* J_xx = dx/dxi */
127: J01 = J01+GNi[0][i]*cy; /* J_xy = dy/dxi */
128: J10 = J10+GNi[1][i]*cx; /* J_yx = dx/deta */
129: J11 = J11+GNi[1][i]*cy; /* J_yy = dy/deta */
130: }
131: J = (J00*J11)-(J01*J10);
133: iJ00 = J11/J;
134: iJ01 = -J01/J;
135: iJ10 = -J10/J;
136: iJ11 = J00/J;
139: for (i = 0; i < NODES_PER_EL; i++) {
140: GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
141: GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
142: }
144: if (det_J != NULL) *det_J = J;
145: }
147: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
148: {
149: *ngp = 4;
150: gp_xi[0][0] = -0.57735026919;gp_xi[0][1] = -0.57735026919;
151: gp_xi[1][0] = -0.57735026919;gp_xi[1][1] = 0.57735026919;
152: gp_xi[2][0] = 0.57735026919;gp_xi[2][1] = 0.57735026919;
153: gp_xi[3][0] = 0.57735026919;gp_xi[3][1] = -0.57735026919;
154: gp_weight[0] = 1.0;
155: gp_weight[1] = 1.0;
156: gp_weight[2] = 1.0;
157: gp_weight[3] = 1.0;
158: }
161: /* procs to the left claim the ghost node as their element */
164: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
165: {
166: PetscInt m,n,p,M,N,P;
167: PetscInt sx,sy,sz;
170: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
171: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
173: if (mxl != NULL) {
174: *mxl = m;
175: if ((sx+m) == M) *mxl = m-1; /* last proc */
176: }
177: if (myl != NULL) {
178: *myl = n;
179: if ((sy+n) == N) *myl = n-1; /* last proc */
180: }
181: if (mzl != NULL) {
182: *mzl = p;
183: if ((sz+p) == P) *mzl = p-1; /* last proc */
184: }
185: return(0);
186: }
190: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
191: {
192: PetscInt si,sj,sk;
195: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
197: *sx = si;
198: if (si) *sx = si+1;
200: *sy = sj;
201: if (sj) *sy = sj+1;
203: if (sk) {
204: *sz = sk;
205: if (sk != 0) *sz = sk+1;
206: }
208: DMDAGetLocalElementSize(da,mx,my,mz);
209: return(0);
210: }
212: /*
213: i,j are the element indices
214: The unknown is a vector quantity.
215: The s[].c is used to indicate the degree of freedom.
216: */
219: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
220: {
222: /* velocity */
223: /* node 0 */
224: s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0; /* Vx0 */
225: s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1; /* Vy0 */
227: /* node 1 */
228: s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0; /* Vx1 */
229: s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1; /* Vy1 */
231: /* node 2 */
232: s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0; /* Vx2 */
233: s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1; /* Vy2 */
235: /* node 3 */
236: s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0; /* Vx3 */
237: s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1; /* Vy3 */
240: /* pressure */
241: s_p[0].i = i;s_p[0].j = j;s_p[0].c = 2; /* P0 */
242: s_p[1].i = i;s_p[1].j = j+1;s_p[1].c = 2; /* P0 */
243: s_p[2].i = i+1;s_p[2].j = j+1;s_p[2].c = 2; /* P1 */
244: s_p[3].i = i+1;s_p[3].j = j;s_p[3].c = 2; /* P1 */
245: return(0);
246: }
250: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
251: {
253: PetscMPIInt rank;
254: PetscInt proc_I,proc_J;
255: PetscInt cpu_x,cpu_y;
256: PetscInt local_mx,local_my;
257: Vec vlx,vly;
258: PetscInt *LX,*LY,i;
259: PetscScalar *_a;
260: Vec V_SEQ;
261: VecScatter ctx;
264: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
266: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
268: proc_J = rank/cpu_x;
269: proc_I = rank-cpu_x*proc_J;
271: PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);
272: PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);
274: DMDAGetLocalElementSize(da,&local_mx,&local_my,NULL);
275: VecCreate(PETSC_COMM_WORLD,&vlx);
276: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
277: VecSetFromOptions(vlx);
279: VecCreate(PETSC_COMM_WORLD,&vly);
280: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
281: VecSetFromOptions(vly);
283: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
284: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
285: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
286: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
290: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
291: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
292: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
293: VecGetArray(V_SEQ,&_a);
294: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
295: VecRestoreArray(V_SEQ,&_a);
296: VecScatterDestroy(&ctx);
297: VecDestroy(&V_SEQ);
299: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
300: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
301: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
302: VecGetArray(V_SEQ,&_a);
303: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
304: VecRestoreArray(V_SEQ,&_a);
305: VecScatterDestroy(&ctx);
306: VecDestroy(&V_SEQ);
310: *_lx = LX;
311: *_ly = LY;
313: VecDestroy(&vlx);
314: VecDestroy(&vly);
315: return(0);
316: }
320: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
321: {
322: DM cda;
323: Vec coords;
324: DMDACoor2d **_coords;
325: PetscInt si,sj,nx,ny,i,j;
326: FILE *fp;
327: char fname[PETSC_MAX_PATH_LEN];
328: PetscMPIInt rank;
332: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
333: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
334: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
335: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
337: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
339: DMGetCoordinateDM(da,&cda);
340: DMGetCoordinatesLocal(da,&coords);
341: DMDAVecGetArray(cda,coords,&_coords);
342: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
343: for (j = sj; j < sj+ny-1; j++) {
344: for (i = si; i < si+nx-1; i++) {
345: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
346: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i].x),PetscRealPart(_coords[j+1][i].y));
347: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i+1].x),PetscRealPart(_coords[j+1][i+1].y));
348: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i+1].x),PetscRealPart(_coords[j][i+1].y));
349: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
350: }
351: }
352: DMDAVecRestoreArray(cda,coords,&_coords);
354: PetscFClose(PETSC_COMM_SELF,fp);
355: return(0);
356: }
360: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
361: {
362: DM cda;
363: Vec coords,local_fields;
364: DMDACoor2d **_coords;
365: FILE *fp;
366: char fname[PETSC_MAX_PATH_LEN];
367: PetscMPIInt rank;
368: PetscInt si,sj,nx,ny,i,j;
369: PetscInt n_dofs,d;
370: PetscScalar *_fields;
374: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
375: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
376: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
377: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
379: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
380: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
381: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
382: for (d = 0; d < n_dofs; d++) {
383: const char *field_name;
384: DMDAGetFieldName(da,d,&field_name);
385: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
386: }
387: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
390: DMGetCoordinateDM(da,&cda);
391: DMGetCoordinatesLocal(da,&coords);
392: DMDAVecGetArray(cda,coords,&_coords);
393: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
395: DMCreateLocalVector(da,&local_fields);
396: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
397: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
398: VecGetArray(local_fields,&_fields);
401: for (j = sj; j < sj+ny; j++) {
402: for (i = si; i < si+nx; i++) {
403: PetscScalar coord_x,coord_y;
404: PetscScalar field_d;
406: coord_x = _coords[j][i].x;
407: coord_y = _coords[j][i].y;
409: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
410: for (d = 0; d < n_dofs; d++) {
411: field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
412: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",PetscRealPart(field_d));
413: }
414: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
415: }
416: }
417: VecRestoreArray(local_fields,&_fields);
418: VecDestroy(&local_fields);
420: DMDAVecRestoreArray(cda,coords,&_coords);
422: PetscFClose(PETSC_COMM_SELF,fp);
423: return(0);
424: }
428: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
429: {
430: DM cda;
431: Vec local_fields;
432: FILE *fp;
433: char fname[PETSC_MAX_PATH_LEN];
434: PetscMPIInt rank;
435: PetscInt si,sj,nx,ny,i,j,p;
436: PetscInt n_dofs,d;
437: GaussPointCoefficients **_coefficients;
438: PetscErrorCode ierr;
441: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
442: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
443: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
444: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
446: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
447: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
448: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
449: for (d = 0; d < n_dofs; d++) {
450: const char *field_name;
451: DMDAGetFieldName(da,d,&field_name);
452: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
453: }
454: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
457: DMGetCoordinateDM(da,&cda);
458: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
460: DMCreateLocalVector(da,&local_fields);
461: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
462: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
463: DMDAVecGetArray(da,local_fields,&_coefficients);
466: for (j = sj; j < sj+ny; j++) {
467: for (i = si; i < si+nx; i++) {
468: PetscScalar coord_x,coord_y;
470: for (p = 0; p < GAUSS_POINTS; p++) {
471: coord_x = _coefficients[j][i].gp_coords[2*p];
472: coord_y = _coefficients[j][i].gp_coords[2*p+1];
474: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
476: 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]));
477: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
478: }
479: }
480: }
481: DMDAVecRestoreArray(da,local_fields,&_coefficients);
482: VecDestroy(&local_fields);
484: PetscFClose(PETSC_COMM_SELF,fp);
485: return(0);
486: }
489: 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)
490: {
491: PetscInt ij;
492: PetscInt r,c,nc;
494: nc = u_NPE*u_dof;
496: r = w_dof*wi+wd;
497: c = u_dof*ui+ud;
499: ij = r*nc+c;
501: return ij;
502: }
504: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
505: {
506: PetscInt ngp;
507: PetscScalar gp_xi[GAUSS_POINTS][2];
508: PetscScalar gp_weight[GAUSS_POINTS];
509: PetscInt p,i,j,k;
510: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
511: PetscScalar J_p,tildeD[3];
512: PetscScalar B[3][U_DOFS*NODES_PER_EL];
515: /* define quadrature rule */
516: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
518: /* evaluate integral */
519: for (p = 0; p < ngp; p++) {
520: ConstructQ12D_GNi(gp_xi[p],GNi_p);
521: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
523: for (i = 0; i < NODES_PER_EL; i++) {
524: PetscScalar d_dx_i = GNx_p[0][i];
525: PetscScalar d_dy_i = GNx_p[1][i];
527: B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
528: B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
529: B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
530: }
533: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
534: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
535: tildeD[2] = gp_weight[p]*J_p*eta[p];
537: /* form Bt tildeD B */
538: /*
539: Ke_ij = Bt_ik . D_kl . B_lj
540: = B_ki . D_kl . B_lj
541: = B_ki . D_kk . B_kj
542: */
543: for (i = 0; i < 8; i++) {
544: for (j = 0; j < 8; j++) {
545: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
546: Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
547: }
548: }
549: }
550: }
551: }
553: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
554: {
555: PetscInt ngp;
556: PetscScalar gp_xi[GAUSS_POINTS][2];
557: PetscScalar gp_weight[GAUSS_POINTS];
558: PetscInt p,i,j,di;
559: PetscScalar Ni_p[NODES_PER_EL];
560: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
561: PetscScalar J_p,fac;
564: /* define quadrature rule */
565: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
567: /* evaluate integral */
568: for (p = 0; p < ngp; p++) {
569: ConstructQ12D_Ni(gp_xi[p],Ni_p);
570: ConstructQ12D_GNi(gp_xi[p],GNi_p);
571: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
572: fac = gp_weight[p]*J_p;
574: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
575: for (di = 0; di < NSD; di++) { /* u dofs */
576: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
577: PetscInt IJ;
578: /* Ke[4*u_idx+j] = Ke[4*u_idx+j] - GNx_p[di][i] * Ni_p[j] * fac; */
579: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);
581: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
582: }
583: }
584: }
585: }
586: }
588: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
589: {
590: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
591: PetscInt i,j;
592: PetscInt nr_g,nc_g;
594: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
595: FormGradientOperatorQ1(Ge,coords);
597: nr_g = U_DOFS*NODES_PER_EL;
598: nc_g = P_DOFS*NODES_PER_EL;
600: for (i = 0; i < nr_g; i++) {
601: for (j = 0; j < nc_g; j++) {
602: De[nr_g*j+i] = Ge[nc_g*i+j];
603: }
604: }
605: }
607: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
608: {
609: PetscInt ngp;
610: PetscScalar gp_xi[GAUSS_POINTS][2];
611: PetscScalar gp_weight[GAUSS_POINTS];
612: PetscInt p,i,j;
613: PetscScalar Ni_p[NODES_PER_EL];
614: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
615: PetscScalar J_p,fac,eta_avg;
618: /* define quadrature rule */
619: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
621: /* evaluate integral */
622: for (p = 0; p < ngp; p++) {
623: ConstructQ12D_Ni(gp_xi[p],Ni_p);
624: ConstructQ12D_GNi(gp_xi[p],GNi_p);
625: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
626: fac = gp_weight[p]*J_p;
628: for (i = 0; i < NODES_PER_EL; i++) {
629: for (j = 0; j < NODES_PER_EL; j++) {
630: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
631: }
632: }
633: }
635: /* scale */
636: eta_avg = 0.0;
637: for (p = 0; p < ngp; p++) eta_avg += eta[p];
638: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
639: fac = 1.0/eta_avg;
640: for (i = 0; i < NODES_PER_EL; i++) {
641: for (j = 0; j < NODES_PER_EL; j++) {
642: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
643: }
644: }
645: }
647: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
648: {
649: PetscInt ngp;
650: PetscScalar gp_xi[GAUSS_POINTS][2];
651: PetscScalar gp_weight[GAUSS_POINTS];
652: PetscInt p,i,j;
653: PetscScalar Ni_p[NODES_PER_EL];
654: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
655: PetscScalar J_p,fac,eta_avg;
658: /* define quadrature rule */
659: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
661: /* evaluate integral */
662: for (p = 0; p < ngp; p++) {
663: ConstructQ12D_Ni(gp_xi[p],Ni_p);
664: ConstructQ12D_GNi(gp_xi[p],GNi_p);
665: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
666: fac = gp_weight[p]*J_p;
668: for (i = 0; i < NODES_PER_EL; i++) {
669: for (j = 0; j < NODES_PER_EL; j++) {
670: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
671: }
672: }
673: }
675: /* scale */
676: eta_avg = 0.0;
677: for (p = 0; p < ngp; p++) eta_avg += eta[p];
678: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
679: fac = 1.0/eta_avg;
680: for (i = 0; i < NODES_PER_EL; i++) {
681: for (j = 0; j < NODES_PER_EL; j++) {
682: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
683: }
684: }
685: }
687: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
688: {
689: PetscInt ngp;
690: PetscScalar gp_xi[GAUSS_POINTS][2];
691: PetscScalar gp_weight[GAUSS_POINTS];
692: PetscInt p,i;
693: PetscScalar Ni_p[NODES_PER_EL];
694: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
695: PetscScalar J_p,fac;
698: /* define quadrature rule */
699: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
701: /* evaluate integral */
702: for (p = 0; p < ngp; p++) {
703: ConstructQ12D_Ni(gp_xi[p],Ni_p);
704: ConstructQ12D_GNi(gp_xi[p],GNi_p);
705: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
706: fac = gp_weight[p]*J_p;
708: for (i = 0; i < NODES_PER_EL; i++) {
709: Fe[NSD*i] += fac*Ni_p[i]*fx[p];
710: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
711: }
712: }
713: }
717: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
718: {
720: /* get coords for the element */
721: el_coords[NSD*0+0] = _coords[ej][ei].x;el_coords[NSD*0+1] = _coords[ej][ei].y;
722: el_coords[NSD*1+0] = _coords[ej+1][ei].x;el_coords[NSD*1+1] = _coords[ej+1][ei].y;
723: el_coords[NSD*2+0] = _coords[ej+1][ei+1].x;el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
724: el_coords[NSD*3+0] = _coords[ej][ei+1].x;el_coords[NSD*3+1] = _coords[ej][ei+1].y;
725: return(0);
726: }
730: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
731: {
732: DM cda;
733: Vec coords;
734: DMDACoor2d **_coords;
735: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
736: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
737: PetscInt sex,sey,mx,my;
738: PetscInt ei,ej;
739: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
740: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
741: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
742: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
743: PetscScalar el_coords[NODES_PER_EL*NSD];
744: Vec local_properties;
745: GaussPointCoefficients **props;
746: PetscScalar *prop_eta;
747: PetscErrorCode ierr;
750: /* setup for coords */
751: DMGetCoordinateDM(stokes_da,&cda);
752: DMGetCoordinatesLocal(stokes_da,&coords);
753: DMDAVecGetArray(cda,coords,&_coords);
755: /* setup for coefficients */
756: DMCreateLocalVector(properties_da,&local_properties);
757: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
758: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
759: DMDAVecGetArray(properties_da,local_properties,&props);
761: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
762: for (ej = sey; ej < sey+my; ej++) {
763: for (ei = sex; ei < sex+mx; ei++) {
764: /* get coords for the element */
765: GetElementCoords(_coords,ei,ej,el_coords);
767: /* get coefficients for the element */
768: prop_eta = props[ej][ei].eta;
770: /* initialise element stiffness matrix */
771: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
772: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
773: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
774: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
776: /* form element stiffness matrix */
777: FormStressOperatorQ1(Ae,el_coords,prop_eta);
778: FormGradientOperatorQ1(Ge,el_coords);
779: FormDivergenceOperatorQ1(De,el_coords);
780: FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);
782: /* insert element matrix into global matrix */
783: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
784: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
785: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
786: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
787: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
788: }
789: }
790: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
791: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
793: DMDAVecRestoreArray(cda,coords,&_coords);
795: DMDAVecRestoreArray(properties_da,local_properties,&props);
796: VecDestroy(&local_properties);
797: return(0);
798: }
802: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
803: {
804: DM cda;
805: Vec coords;
806: DMDACoor2d **_coords;
807: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
808: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
809: PetscInt sex,sey,mx,my;
810: PetscInt ei,ej;
811: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
812: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
813: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
814: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
815: PetscScalar el_coords[NODES_PER_EL*NSD];
816: Vec local_properties;
817: GaussPointCoefficients **props;
818: PetscScalar *prop_eta;
819: PetscErrorCode ierr;
822: /* setup for coords */
823: DMGetCoordinateDM(stokes_da,&cda);
824: DMGetCoordinatesLocal(stokes_da,&coords);
825: DMDAVecGetArray(cda,coords,&_coords);
827: /* setup for coefficients */
828: DMCreateLocalVector(properties_da,&local_properties);
829: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
830: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
831: DMDAVecGetArray(properties_da,local_properties,&props);
833: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
834: for (ej = sey; ej < sey+my; ej++) {
835: for (ei = sex; ei < sex+mx; ei++) {
836: /* get coords for the element */
837: GetElementCoords(_coords,ei,ej,el_coords);
839: /* get coefficients for the element */
840: prop_eta = props[ej][ei].eta;
842: /* initialise element stiffness matrix */
843: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
844: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
845: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
846: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
849: /* form element stiffness matrix */
850: FormStressOperatorQ1(Ae,el_coords,prop_eta);
851: FormGradientOperatorQ1(Ge,el_coords);
852: /* FormDivergenceOperatorQ1(De, el_coords); */
853: FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);
855: /* insert element matrix into global matrix */
856: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
857: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
858: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
859: /* MatSetValuesStencil(A, NODES_PER_EL*P_DOFS,p_eqn, NODES_PER_EL*U_DOFS,u_eqn, De, ADD_VALUES); */
860: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
861: }
862: }
863: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
864: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
866: DMDAVecRestoreArray(cda,coords,&_coords);
868: DMDAVecRestoreArray(properties_da,local_properties,&props);
869: VecDestroy(&local_properties);
870: return(0);
871: }
875: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
876: {
877: PetscInt n;
880: for (n = 0; n < 4; n++) {
881: 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];
882: 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];
883: 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];
884: }
885: return(0);
886: }
890: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
891: {
892: DM cda;
893: Vec coords;
894: DMDACoor2d **_coords;
895: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
896: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
897: PetscInt sex,sey,mx,my;
898: PetscInt ei,ej;
899: PetscScalar Fe[NODES_PER_EL*U_DOFS];
900: PetscScalar He[NODES_PER_EL*P_DOFS];
901: PetscScalar el_coords[NODES_PER_EL*NSD];
902: Vec local_properties;
903: GaussPointCoefficients **props;
904: PetscScalar *prop_fx,*prop_fy;
905: Vec local_F;
906: StokesDOF **ff;
907: PetscErrorCode ierr;
910: /* setup for coords */
911: DMGetCoordinateDM(stokes_da,&cda);
912: DMGetCoordinatesLocal(stokes_da,&coords);
913: DMDAVecGetArray(cda,coords,&_coords);
915: /* setup for coefficients */
916: DMGetLocalVector(properties_da,&local_properties);
917: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
918: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
919: DMDAVecGetArray(properties_da,local_properties,&props);
921: /* get acces to the vector */
922: DMGetLocalVector(stokes_da,&local_F);
923: VecZeroEntries(local_F);
924: DMDAVecGetArray(stokes_da,local_F,&ff);
927: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
928: for (ej = sey; ej < sey+my; ej++) {
929: for (ei = sex; ei < sex+mx; ei++) {
930: /* get coords for the element */
931: GetElementCoords(_coords,ei,ej,el_coords);
933: /* get coefficients for the element */
934: prop_fx = props[ej][ei].fx;
935: prop_fy = props[ej][ei].fy;
937: /* initialise element stiffness matrix */
938: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
939: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
942: /* form element stiffness matrix */
943: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
945: /* insert element matrix into global matrix */
946: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
948: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
949: }
950: }
952: DMDAVecRestoreArray(stokes_da,local_F,&ff);
953: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
954: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
955: DMRestoreLocalVector(stokes_da,&local_F);
958: DMDAVecRestoreArray(cda,coords,&_coords);
960: DMDAVecRestoreArray(properties_da,local_properties,&props);
961: DMRestoreLocalVector(properties_da,&local_properties);
962: return(0);
963: }
967: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM *_da,Vec *_X)
968: {
969: DM da,cda;
970: Vec X,local_X;
971: StokesDOF **_stokes;
972: Vec coords;
973: DMDACoor2d **_coords;
974: PetscInt si,sj,ei,ej,i,j;
978: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
979: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,&da);
980: DMDASetFieldName(da,0,"anlytic_Vx");
981: DMDASetFieldName(da,1,"anlytic_Vy");
982: DMDASetFieldName(da,2,"analytic_P");
985: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.,0.);
988: DMGetCoordinatesLocal(da,&coords);
989: DMGetCoordinateDM(da,&cda);
990: DMDAVecGetArray(cda,coords,&_coords);
992: DMCreateGlobalVector(da,&X);
993: DMCreateLocalVector(da,&local_X);
994: DMDAVecGetArray(da,local_X,&_stokes);
996: DMDAGetGhostCorners(da,&si,&sj,0,&ei,&ej,0);
997: for (j = sj; j < sj+ej; j++) {
998: for (i = si; i < si+ei; i++) {
999: double pos[2],pressure,vel[2],total_stress[3],strain_rate[3];
1001: pos[0] = PetscRealPart(_coords[j][i].x);
1002: pos[1] = PetscRealPart(_coords[j][i].y);
1004: evaluate_solCx(pos,eta0,eta1,xc,nz,vel,&pressure,total_stress,strain_rate);
1006: _stokes[j][i].u_dof = vel[0];
1007: _stokes[j][i].v_dof = vel[1];
1008: _stokes[j][i].p_dof = pressure;
1009: }
1010: }
1011: DMDAVecRestoreArray(da,local_X,&_stokes);
1012: DMDAVecRestoreArray(cda,coords,&_coords);
1014: DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1015: DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);
1017: VecDestroy(&local_X);
1019: *_da = da;
1020: *_X = X;
1021: return(0);
1022: }
1026: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
1027: {
1029: /* get the nodal fields */
1030: 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;
1031: 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;
1032: 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;
1033: 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;
1034: return(0);
1035: }
1039: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
1040: {
1041: DM cda;
1042: Vec coords,X_analytic_local,X_local;
1043: DMDACoor2d **_coords;
1044: PetscInt sex,sey,mx,my;
1045: PetscInt ei,ej;
1046: PetscScalar el_coords[NODES_PER_EL*NSD];
1047: StokesDOF **stokes_analytic,**stokes;
1048: StokesDOF stokes_analytic_e[4],stokes_e[4];
1050: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1051: PetscScalar Ni_p[NODES_PER_EL];
1052: PetscInt ngp;
1053: PetscScalar gp_xi[GAUSS_POINTS][2];
1054: PetscScalar gp_weight[GAUSS_POINTS];
1055: PetscInt p,i;
1056: PetscScalar J_p,fac;
1057: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1058: PetscInt M;
1059: PetscReal xymin[2],xymax[2];
1063: /* define quadrature rule */
1064: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1066: /* setup for coords */
1067: DMGetCoordinateDM(stokes_da,&cda);
1068: DMGetCoordinatesLocal(stokes_da,&coords);
1069: DMDAVecGetArray(cda,coords,&_coords);
1071: /* setup for analytic */
1072: DMCreateLocalVector(stokes_da,&X_analytic_local);
1073: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1074: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1075: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1077: /* setup for solution */
1078: DMCreateLocalVector(stokes_da,&X_local);
1079: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1080: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1081: DMDAVecGetArray(stokes_da,X_local,&stokes);
1083: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1084: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1086: h = (xymax[0]-xymin[0])/((double)M);
1088: tp_L2 = tu_L2 = tu_H1 = 0.0;
1090: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
1091: for (ej = sey; ej < sey+my; ej++) {
1092: for (ei = sex; ei < sex+mx; ei++) {
1093: /* get coords for the element */
1094: GetElementCoords(_coords,ei,ej,el_coords);
1095: StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1096: StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);
1098: /* evaluate integral */
1099: p_e_L2 = 0.0;
1100: u_e_L2 = 0.0;
1101: u_e_H1 = 0.0;
1102: for (p = 0; p < ngp; p++) {
1103: ConstructQ12D_Ni(gp_xi[p],Ni_p);
1104: ConstructQ12D_GNi(gp_xi[p],GNi_p);
1105: ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1106: fac = gp_weight[p]*J_p;
1108: for (i = 0; i < NODES_PER_EL; i++) {
1109: PetscScalar u_error,v_error;
1111: 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);
1113: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1114: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1115: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);
1117: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1118: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error /* du/dy */
1119: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1120: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error); /* dv/dy */
1121: }
1122: }
1124: tp_L2 += p_e_L2;
1125: tu_L2 += u_e_L2;
1126: tu_H1 += u_e_H1;
1127: }
1128: }
1129: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1130: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1131: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1132: p_L2 = PetscSqrtScalar(p_L2);
1133: u_L2 = PetscSqrtScalar(u_L2);
1134: u_H1 = PetscSqrtScalar(u_H1);
1136: PetscPrintf(PETSC_COMM_WORLD,"%1.4e %1.4e %1.4e %1.4e \n",PetscRealPart(h),PetscRealPart(p_L2),PetscRealPart(u_L2),PetscRealPart(u_H1));
1139: DMDAVecRestoreArray(cda,coords,&_coords);
1141: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1142: VecDestroy(&X_analytic_local);
1143: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1144: VecDestroy(&X_local);
1145: return(0);
1146: }
1150: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1151: {
1152: DM da_Stokes,da_prop;
1153: PetscInt u_dof,p_dof,dof,stencil_width;
1154: Mat A,B;
1155: PetscInt mxl,myl;
1156: DM prop_cda,vel_cda;
1157: Vec prop_coords,vel_coords;
1158: PetscInt si,sj,nx,ny,i,j,p;
1159: Vec f,X;
1160: PetscInt prop_dof,prop_stencil_width;
1161: Vec properties,l_properties;
1162: PetscReal dx,dy;
1163: PetscInt M,N;
1164: DMDACoor2d **_prop_coords,**_vel_coords;
1165: GaussPointCoefficients **element_props;
1166: PetscInt its;
1167: KSP ksp_S;
1168: PetscInt coefficient_structure = 0;
1169: PetscInt cpu_x,cpu_y,*lx = NULL,*ly = NULL;
1170: PetscBool use_gp_coords = PETSC_FALSE,set;
1171: char filename[PETSC_MAX_PATH_LEN];
1172: PetscErrorCode ierr;
1175: /* Generate the da for velocity and pressure */
1176: /*
1177: We use Q1 elements for the temperature.
1178: FEM has a 9-point stencil (BOX) or connectivity pattern
1179: Num nodes in each direction is mx+1, my+1
1180: */
1181: u_dof = U_DOFS; /* Vx, Vy - velocities */
1182: p_dof = P_DOFS; /* p - pressure */
1183: dof = u_dof+p_dof;
1184: stencil_width = 1;
1185: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1186: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&da_Stokes);
1187: DMDASetFieldName(da_Stokes,0,"Vx");
1188: DMDASetFieldName(da_Stokes,1,"Vy");
1189: DMDASetFieldName(da_Stokes,2,"P");
1191: /* unit box [0,1] x [0,1] */
1192: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.,0.);
1195: /* Generate element properties, we will assume all material properties are constant over the element */
1196: /* local number of elements */
1197: DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,NULL);
1199: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1200: DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1201: DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);
1203: prop_dof = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1204: prop_stencil_width = 0;
1205: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1206: mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
1207: PetscFree(lx);
1208: PetscFree(ly);
1210: /* define centroid positions */
1211: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1212: dx = 1.0/((PetscReal)(M));
1213: dy = 1.0/((PetscReal)(N));
1215: 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);
1217: /* define coefficients */
1218: PetscOptionsGetInt(NULL,"-c_str",&coefficient_structure,NULL);
1219: /* PetscPrintf(PETSC_COMM_WORLD, "Using coeficient structure %D \n", coefficient_structure); */
1221: DMCreateGlobalVector(da_prop,&properties);
1222: DMCreateLocalVector(da_prop,&l_properties);
1223: DMDAVecGetArray(da_prop,l_properties,&element_props);
1225: DMGetCoordinateDM(da_prop,&prop_cda);
1226: DMGetCoordinatesLocal(da_prop,&prop_coords);
1227: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
1229: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
1231: DMGetCoordinateDM(da_Stokes,&vel_cda);
1232: DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1233: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1236: /* interpolate the coordinates */
1237: for (j = sj; j < sj+ny; j++) {
1238: for (i = si; i < si+nx; i++) {
1239: PetscInt ngp;
1240: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1241: PetscScalar el_coords[8];
1243: GetElementCoords(_vel_coords,i,j,el_coords);
1244: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1246: for (p = 0; p < GAUSS_POINTS; p++) {
1247: PetscScalar gp_x,gp_y;
1248: PetscInt n;
1249: PetscScalar xi_p[2],Ni_p[4];
1251: xi_p[0] = gp_xi[p][0];
1252: xi_p[1] = gp_xi[p][1];
1253: ConstructQ12D_Ni(xi_p,Ni_p);
1255: gp_x = 0.0;
1256: gp_y = 0.0;
1257: for (n = 0; n < NODES_PER_EL; n++) {
1258: gp_x = gp_x+Ni_p[n]*el_coords[2*n];
1259: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1260: }
1261: element_props[j][i].gp_coords[2*p] = gp_x;
1262: element_props[j][i].gp_coords[2*p+1] = gp_y;
1263: }
1264: }
1265: }
1267: /* define the coefficients */
1268: PetscOptionsGetBool(NULL,"-use_gp_coords",&use_gp_coords,0);
1270: for (j = sj; j < sj+ny; j++) {
1271: for (i = si; i < si+nx; i++) {
1272: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1273: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1274: PetscReal coord_x,coord_y;
1276: if (coefficient_structure == 0) {
1277: PetscReal opts_eta0,opts_eta1,opts_xc;
1278: PetscInt opts_nz;
1280: opts_eta0 = 1.0;
1281: opts_eta1 = 1.0;
1282: opts_xc = 0.5;
1283: opts_nz = 1;
1285: PetscOptionsGetReal(NULL,"-solcx_eta0",&opts_eta0,0);
1286: PetscOptionsGetReal(NULL,"-solcx_eta1",&opts_eta1,0);
1287: PetscOptionsGetReal(NULL,"-solcx_xc",&opts_xc,0);
1288: PetscOptionsGetInt(NULL,"-solcx_nz",&opts_nz,0);
1290: for (p = 0; p < GAUSS_POINTS; p++) {
1291: coord_x = centroid_x;
1292: coord_y = centroid_y;
1293: if (use_gp_coords) {
1294: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1295: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1296: }
1299: element_props[j][i].eta[p] = opts_eta0;
1300: if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;
1302: element_props[j][i].fx[p] = 0.0;
1303: element_props[j][i].fy[p] = sin((double)opts_nz*PETSC_PI*coord_y)*cos(1.0*PETSC_PI*coord_x);
1304: }
1305: } else if (coefficient_structure == 1) { /* square sinker */
1306: PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;
1308: opts_eta0 = 1.0;
1309: opts_eta1 = 1.0;
1310: opts_dx = 0.50;
1311: opts_dy = 0.50;
1313: PetscOptionsGetReal(NULL,"-sinker_eta0",&opts_eta0,0);
1314: PetscOptionsGetReal(NULL,"-sinker_eta1",&opts_eta1,0);
1315: PetscOptionsGetReal(NULL,"-sinker_dx",&opts_dx,0);
1316: PetscOptionsGetReal(NULL,"-sinker_dy",&opts_dy,0);
1319: for (p = 0; p < GAUSS_POINTS; p++) {
1320: coord_x = centroid_x;
1321: coord_y = centroid_y;
1322: if (use_gp_coords) {
1323: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1324: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1325: }
1327: element_props[j][i].eta[p] = opts_eta0;
1328: element_props[j][i].fx[p] = 0.0;
1329: element_props[j][i].fy[p] = 0.0;
1331: if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1332: if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1333: element_props[j][i].eta[p] = opts_eta1;
1334: element_props[j][i].fx[p] = 0.0;
1335: element_props[j][i].fy[p] = -1.0;
1336: }
1337: }
1338: }
1339: } else if (coefficient_structure == 2) { /* circular sinker */
1340: PetscReal opts_eta0,opts_eta1,opts_r,radius2;
1342: opts_eta0 = 1.0;
1343: opts_eta1 = 1.0;
1344: opts_r = 0.25;
1346: PetscOptionsGetReal(NULL,"-sinker_eta0",&opts_eta0,0);
1347: PetscOptionsGetReal(NULL,"-sinker_eta1",&opts_eta1,0);
1348: PetscOptionsGetReal(NULL,"-sinker_r",&opts_r,0);
1350: for (p = 0; p < GAUSS_POINTS; p++) {
1351: coord_x = centroid_x;
1352: coord_y = centroid_y;
1353: if (use_gp_coords) {
1354: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1355: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1356: }
1358: element_props[j][i].eta[p] = opts_eta0;
1359: element_props[j][i].fx[p] = 0.0;
1360: element_props[j][i].fy[p] = 0.0;
1362: radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1363: if (radius2 < opts_r*opts_r) {
1364: element_props[j][i].eta[p] = opts_eta1;
1365: element_props[j][i].fx[p] = 0.0;
1366: element_props[j][i].fy[p] = -1.0;
1367: }
1368: }
1369: } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1370: PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;
1372: opts_eta0 = 1.0;
1373: opts_eta1 = 1.0;
1374: opts_r = 0.25;
1375: opts_c0x = 0.35; /* circle center */
1376: opts_c0y = 0.35;
1377: opts_s0x = 0.7; /* square center */
1378: opts_s0y = 0.7;
1379: opts_dx = 0.25;
1380: opts_dy = 0.25;
1381: opts_phi = 25;
1383: PetscOptionsGetReal(NULL,"-sinker_eta0",&opts_eta0,0);
1384: PetscOptionsGetReal(NULL,"-sinker_eta1",&opts_eta1,0);
1385: PetscOptionsGetReal(NULL,"-sinker_r",&opts_r,0);
1386: PetscOptionsGetReal(NULL,"-sinker_c0x",&opts_c0x,0);
1387: PetscOptionsGetReal(NULL,"-sinker_c0y",&opts_c0y,0);
1388: PetscOptionsGetReal(NULL,"-sinker_s0x",&opts_s0x,0);
1389: PetscOptionsGetReal(NULL,"-sinker_s0y",&opts_s0y,0);
1390: PetscOptionsGetReal(NULL,"-sinker_dx",&opts_dx,0);
1391: PetscOptionsGetReal(NULL,"-sinker_dy",&opts_dy,0);
1392: PetscOptionsGetReal(NULL,"-sinker_phi",&opts_phi,0);
1393: opts_phi *= PETSC_PI / 180;
1395: for (p = 0; p < GAUSS_POINTS; p++) {
1396: coord_x = centroid_x;
1397: coord_y = centroid_y;
1398: if (use_gp_coords) {
1399: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1400: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1401: }
1403: element_props[j][i].eta[p] = opts_eta0;
1404: element_props[j][i].fx[p] = 0.0;
1405: element_props[j][i].fy[p] = -0.2;
1407: radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1408: if (radius2 < opts_r*opts_r
1409: || (PetscAbs(+(coord_x - opts_s0x)*cos(opts_phi) + (coord_y - opts_s0y)*sin(opts_phi)) < opts_dx/2
1410: && PetscAbs(-(coord_x - opts_s0x)*sin(opts_phi) + (coord_y - opts_s0y)*cos(opts_phi)) < opts_dy/2)) {
1411: element_props[j][i].eta[p] = opts_eta1;
1412: element_props[j][i].fx[p] = 0.0;
1413: element_props[j][i].fy[p] = -1.0;
1414: }
1415: }
1416: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1417: }
1418: }
1419: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
1421: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1423: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1424: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1425: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
1428: DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1429: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for Stokes eqn.","properties");
1432: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1433: DMCreateMatrix(da_Stokes,MATAIJ,&A);
1434: DMCreateMatrix(da_Stokes,MATAIJ,&B);
1435: DMCreateGlobalVector(da_Stokes,&f);
1436: DMCreateGlobalVector(da_Stokes,&X);
1438: /* assemble A11 */
1439: MatZeroEntries(A);
1440: MatZeroEntries(B);
1441: VecZeroEntries(f);
1443: AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1444: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1445: /* build force vector */
1446: AssembleF_Stokes(f,da_Stokes,da_prop,properties);
1448: DMDABCApplyFreeSlip(da_Stokes,A,f);
1449: DMDABCApplyFreeSlip(da_Stokes,B,NULL);
1451: /* SOLVE */
1452: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1453: KSPSetOptionsPrefix(ksp_S,"stokes_");
1454: KSPSetOperators(ksp_S,A,B,SAME_NONZERO_PATTERN);
1455: KSPSetDM(ksp_S,da_Stokes);
1456: KSPSetDMActive(ksp_S,PETSC_FALSE);
1457: KSPSetFromOptions(ksp_S);
1458: {
1459: PC pc;
1460: const PetscInt ufields[] = {0,1},pfields[1] = {2};
1461: KSPGetPC(ksp_S,&pc);
1462: PCFieldSplitSetBlockSize(pc,3);
1463: PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1464: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1465: }
1467: KSPSolve(ksp_S,f,X);
1469: PetscOptionsGetString(NULL,"-o",filename,sizeof(filename),&set);
1470: if (set) {
1471: char *ext;
1472: PetscViewer viewer;
1473: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1474: PetscStrrchr(filename,'.',&ext);
1475: if (!strcmp("vts",ext)) {
1476: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1477: } else {
1478: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1479: }
1480: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1481: PetscViewerFileSetName(viewer,filename);
1482: VecView(X,viewer);
1483: PetscViewerDestroy(&viewer);
1484: }
1485: DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
1487: KSPGetIterationNumber(ksp_S,&its);
1489: /*
1490: {
1491: Vec x;
1492: PetscInt L;
1493: VecDuplicate(f,&x);
1494: MatMult(A,X, x);
1495: VecAXPY(x, -1.0, f);
1496: VecNorm(x, NORM_2, &nrm);
1497: PetscPrintf(PETSC_COMM_WORLD, "Its. %1.4d, norm |AX-f| = %1.5e \n", its, nrm);
1498: VecDestroy(&x);
1500: VecNorm(X, NORM_2, &nrm);
1501: VecGetSize(X, &L);
1502: PetscPrintf(PETSC_COMM_WORLD, " norm |X|/sqrt{N} = %1.5e \n", nrm/sqrt((PetscScalar)L));
1503: }
1504: */
1506: if (coefficient_structure == 0) {
1507: PetscReal opts_eta0,opts_eta1,opts_xc;
1508: PetscInt opts_nz,N;
1509: DM da_Stokes_analytic;
1510: Vec X_analytic;
1511: PetscReal nrm1[3],nrm2[3],nrmI[3];
1513: opts_eta0 = 1.0;
1514: opts_eta1 = 1.0;
1515: opts_xc = 0.5;
1516: opts_nz = 1;
1518: PetscOptionsGetReal(NULL,"-solcx_eta0",&opts_eta0,0);
1519: PetscOptionsGetReal(NULL,"-solcx_eta1",&opts_eta1,0);
1520: PetscOptionsGetReal(NULL,"-solcx_xc",&opts_xc,0);
1521: PetscOptionsGetInt(NULL,"-solcx_nz",&opts_nz,0);
1524: DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1525: DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
1527: DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);
1530: VecAXPY(X_analytic,-1.0,X);
1531: VecGetSize(X_analytic,&N);
1532: N = N/3;
1534: VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1535: VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1536: VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);
1538: VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1539: VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1540: VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);
1542: VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1543: VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1544: VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);
1546: DMDestroy(&da_Stokes_analytic);
1547: VecDestroy(&X_analytic);
1548: }
1551: KSPDestroy(&ksp_S);
1552: VecDestroy(&X);
1553: VecDestroy(&f);
1554: MatDestroy(&A);
1555: MatDestroy(&B);
1557: DMDestroy(&da_Stokes);
1558: DMDestroy(&da_prop);
1560: VecDestroy(&properties);
1561: VecDestroy(&l_properties);
1562: return(0);
1563: }
1567: int main(int argc,char **args)
1568: {
1570: PetscInt mx,my;
1572: PetscInitialize(&argc,&args,(char*)0,help);
1574: mx = my = 10;
1575: PetscOptionsGetInt(NULL,"-mx",&mx,NULL);
1576: PetscOptionsGetInt(NULL,"-my",&my,NULL);
1578: solve_stokes_2d_coupled(mx,my);
1580: PetscFinalize();
1581: return 0;
1582: }
1584: /* -------------------------- helpers for boundary conditions -------------------------------- */
1588: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1589: {
1590: DM cda;
1591: Vec coords;
1592: PetscInt si,sj,nx,ny,i,j;
1593: PetscInt M,N;
1594: DMDACoor2d **_coords;
1595: PetscInt *g_idx;
1596: PetscInt *bc_global_ids;
1597: PetscScalar *bc_vals;
1598: PetscInt nbcs;
1599: PetscInt n_dofs;
1603: /* enforce bc's */
1604: DMDAGetGlobalIndices(da,NULL,&g_idx);
1606: DMGetCoordinateDM(da,&cda);
1607: DMGetCoordinatesLocal(da,&coords);
1608: DMDAVecGetArray(cda,coords,&_coords);
1609: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1610: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1612: /* --- */
1614: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1615: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1617: /* init the entries to -1 so VecSetValues will ignore them */
1618: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1620: i = nx-1;
1621: for (j = 0; j < ny; j++) {
1622: PetscInt local_id;
1624: local_id = i+j*nx;
1626: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1628: bc_vals[j] = bc_val;
1629: }
1630: nbcs = 0;
1631: if ((si+nx) == (M)) nbcs = ny;
1633: if (b != NULL) {
1634: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1635: VecAssemblyBegin(b);
1636: VecAssemblyEnd(b);
1637: }
1638: if (A != NULL) {
1639: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1640: }
1643: PetscFree(bc_vals);
1644: PetscFree(bc_global_ids);
1646: DMDAVecRestoreArray(cda,coords,&_coords);
1647: return(0);
1648: }
1652: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1653: {
1654: DM cda;
1655: Vec coords;
1656: PetscInt si,sj,nx,ny,i,j;
1657: PetscInt M,N;
1658: DMDACoor2d **_coords;
1659: PetscInt *g_idx;
1660: PetscInt *bc_global_ids;
1661: PetscScalar *bc_vals;
1662: PetscInt nbcs;
1663: PetscInt n_dofs;
1667: /* enforce bc's */
1668: DMDAGetGlobalIndices(da,NULL,&g_idx);
1670: DMGetCoordinateDM(da,&cda);
1671: DMGetCoordinatesLocal(da,&coords);
1672: DMDAVecGetArray(cda,coords,&_coords);
1673: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1674: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1676: /* --- */
1678: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1679: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1681: /* init the entries to -1 so VecSetValues will ignore them */
1682: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1684: i = 0;
1685: for (j = 0; j < ny; j++) {
1686: PetscInt local_id;
1688: local_id = i+j*nx;
1690: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1692: bc_vals[j] = bc_val;
1693: }
1694: nbcs = 0;
1695: if (si == 0) nbcs = ny;
1697: if (b != NULL) {
1698: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1699: VecAssemblyBegin(b);
1700: VecAssemblyEnd(b);
1701: }
1702: if (A != NULL) {
1703: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1704: }
1707: PetscFree(bc_vals);
1708: PetscFree(bc_global_ids);
1710: DMDAVecRestoreArray(cda,coords,&_coords);
1711: return(0);
1712: }
1716: static PetscErrorCode BCApply_NORTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1717: {
1718: DM cda;
1719: Vec coords;
1720: PetscInt si,sj,nx,ny,i,j;
1721: PetscInt M,N;
1722: DMDACoor2d **_coords;
1723: PetscInt *g_idx;
1724: PetscInt *bc_global_ids;
1725: PetscScalar *bc_vals;
1726: PetscInt nbcs;
1727: PetscInt n_dofs;
1731: /* enforce bc's */
1732: DMDAGetGlobalIndices(da,NULL,&g_idx);
1734: DMGetCoordinateDM(da,&cda);
1735: DMGetCoordinatesLocal(da,&coords);
1736: DMDAVecGetArray(cda,coords,&_coords);
1737: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1738: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1740: /* --- */
1742: PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1743: PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);
1745: /* init the entries to -1 so VecSetValues will ignore them */
1746: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1748: j = ny-1;
1749: for (i = 0; i < nx; i++) {
1750: PetscInt local_id;
1752: local_id = i+j*nx;
1754: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1756: bc_vals[i] = bc_val;
1757: }
1758: nbcs = 0;
1759: if ((sj+ny) == (N)) nbcs = nx;
1761: if (b != NULL) {
1762: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1763: VecAssemblyBegin(b);
1764: VecAssemblyEnd(b);
1765: }
1766: if (A != NULL) {
1767: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1768: }
1771: PetscFree(bc_vals);
1772: PetscFree(bc_global_ids);
1774: DMDAVecRestoreArray(cda,coords,&_coords);
1775: return(0);
1776: }
1780: static PetscErrorCode BCApply_SOUTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1781: {
1782: DM cda;
1783: Vec coords;
1784: PetscInt si,sj,nx,ny,i,j;
1785: PetscInt M,N;
1786: DMDACoor2d **_coords;
1787: PetscInt *g_idx;
1788: PetscInt *bc_global_ids;
1789: PetscScalar *bc_vals;
1790: PetscInt nbcs;
1791: PetscInt n_dofs;
1795: /* enforce bc's */
1796: DMDAGetGlobalIndices(da,NULL,&g_idx);
1798: DMGetCoordinateDM(da,&cda);
1799: DMGetCoordinatesLocal(da,&coords);
1800: DMDAVecGetArray(cda,coords,&_coords);
1801: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1802: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1804: /* --- */
1806: PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1807: PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);
1809: /* init the entries to -1 so VecSetValues will ignore them */
1810: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1812: j = 0;
1813: for (i = 0; i < nx; i++) {
1814: PetscInt local_id;
1816: local_id = i+j*nx;
1818: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1820: bc_vals[i] = bc_val;
1821: }
1822: nbcs = 0;
1823: if (sj == 0) nbcs = nx;
1825: if (b != NULL) {
1826: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1827: VecAssemblyBegin(b);
1828: VecAssemblyEnd(b);
1829: }
1830: if (A != NULL) {
1831: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1832: }
1835: PetscFree(bc_vals);
1836: PetscFree(bc_global_ids);
1838: DMDAVecRestoreArray(cda,coords,&_coords);
1839: return(0);
1840: }
1844: /*
1845: Free slip sides.
1846: */
1849: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1850: {
1854: BCApply_NORTH(da_Stokes,1,0.0,A,f);
1855: BCApply_EAST(da_Stokes,0,0.0,A,f);
1856: BCApply_SOUTH(da_Stokes,1,0.0,A,f);
1857: BCApply_WEST(da_Stokes,0,0.0,A,f);
1858: return(0);
1859: }