Actual source code: ex43.c
petsc-3.7.5 2017-01-01
1: static char help[] = "Solves the incompressible, variable viscosity Stokes equation in 2d on the unit domain \n\
2: using Q1Q1 elements, stabilized with Bochev's polynomial projection method. \n\
3: The models defined utilise free slip boundary conditions on all sides. \n\
4: Options: \n"
5: "\
6: -mx : Number of elements in the x-direction \n\
7: -my : Number of elements in the y-direction \n\
8: -o : Specify output filename for solution (will be petsc binary format or paraview format if the extension is .vts) \n\
9: -gnuplot : Output Gauss point coordinates, coefficients and u,p solution in gnuplot format \n\
10: -c_str : Indicates the structure of the coefficients to use \n"
11: "\
12: -c_str 0 => Coefficient definition for an analytic solution with a vertical jump in viscosity at x = xc \n\
13: This problem is driven by the forcing function f(x,y) = (0, sin(nz pi y)cos(pi x) \n\
14: Parameters: \n\
15: -solcx_eta0 : Viscosity to the left of the interface \n\
16: -solcx_eta1 : Viscosity to the right of the interface \n\
17: -solcx_xc : Location of the interface \n\
18: -solcx_nz : Wavenumber in the y direction \n"
19: "\
20: -c_str 1 => Coefficient definition for a dense rectangular blob located at the center of the domain \n\
21: Parameters: \n\
22: -sinker_eta0 : Viscosity of the background fluid \n\
23: -sinker_eta1 : Viscosity of the blob \n\
24: -sinker_dx : Width of the blob \n\
25: -sinker_dy : Height of the blob \n"
26: "\
27: -c_str 2 => Coefficient definition for a dense circular blob located at the center of the domain \n\
28: Parameters: \n\
29: -sinker_eta0 : Viscosity of the background fluid \n\
30: -sinker_eta1 : Viscosity of the blob \n\
31: -sinker_r : Radius of the blob \n"
32: "\
33: -c_str 3 => Coefficient definition for a dense circular and rectangular inclusion (located at the center of the domain) \n\
34: -sinker_eta0 : Viscosity of the background fluid \n\
35: -sinker_eta1 : Viscosity of the two inclusions \n\
36: -sinker_r : Radius of the circular inclusion \n\
37: -sinker_c0x : Origin (x-coord) of the circular inclusion \n\
38: -sinker_c0y : Origin (y-coord) of the circular inclusion \n\
39: -sinker_dx : Width of the rectangular inclusion \n\
40: -sinker_dy : Height of the rectangular inclusion \n\
41: -sinker_phi : Rotation angle of the rectangular inclusion \n\
42: -use_gp_coords : Evaluate the viscosity and force term at the global coordinates of each quadrature point \n\
43: By default, the viscosity and force term are evaulated at the element center and applied as a constant over the entire element \n";
45: /* Contributed by Dave May */
47: #include <petscksp.h>
48: #include <petscdm.h>
49: #include <petscdmda.h>
51: /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
52: #include "ex43-solcx.h"
54: static PetscErrorCode DMDABCApplyFreeSlip(DM,Mat,Vec);
57: #define NSD 2 /* number of spatial dimensions */
58: #define NODES_PER_EL 4 /* nodes per element */
59: #define U_DOFS 2 /* degrees of freedom per velocity node */
60: #define P_DOFS 1 /* degrees of freedom per pressure node */
61: #define GAUSS_POINTS 4
63: /* Gauss point based evaluation 8+4+4+4 = 20 */
64: typedef struct {
65: PetscScalar gp_coords[2*GAUSS_POINTS];
66: PetscScalar eta[GAUSS_POINTS];
67: PetscScalar fx[GAUSS_POINTS];
68: PetscScalar fy[GAUSS_POINTS];
69: } GaussPointCoefficients;
71: typedef struct {
72: PetscScalar u_dof;
73: PetscScalar v_dof;
74: PetscScalar p_dof;
75: } StokesDOF;
78: /*
79: Element: Local basis function ordering
80: 1-----2
81: | |
82: | |
83: 0-----3
84: */
85: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
86: {
87: PetscScalar xi = _xi[0];
88: PetscScalar eta = _xi[1];
90: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
91: Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
92: Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
93: Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
94: }
96: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
97: {
98: PetscScalar xi = _xi[0];
99: PetscScalar eta = _xi[1];
101: GNi[0][0] = -0.25*(1.0-eta);
102: GNi[0][1] = -0.25*(1.0+eta);
103: GNi[0][2] = 0.25*(1.0+eta);
104: GNi[0][3] = 0.25*(1.0-eta);
106: GNi[1][0] = -0.25*(1.0-xi);
107: GNi[1][1] = 0.25*(1.0-xi);
108: GNi[1][2] = 0.25*(1.0+xi);
109: GNi[1][3] = -0.25*(1.0+xi);
110: }
112: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
113: {
114: PetscScalar J00,J01,J10,J11,J;
115: PetscScalar iJ00,iJ01,iJ10,iJ11;
116: PetscInt i;
118: J00 = J01 = J10 = J11 = 0.0;
119: for (i = 0; i < NODES_PER_EL; i++) {
120: PetscScalar cx = coords[2*i];
121: PetscScalar cy = coords[2*i+1];
123: J00 = J00+GNi[0][i]*cx; /* J_xx = dx/dxi */
124: J01 = J01+GNi[0][i]*cy; /* J_xy = dy/dxi */
125: J10 = J10+GNi[1][i]*cx; /* J_yx = dx/deta */
126: J11 = J11+GNi[1][i]*cy; /* J_yy = dy/deta */
127: }
128: J = (J00*J11)-(J01*J10);
130: iJ00 = J11/J;
131: iJ01 = -J01/J;
132: iJ10 = -J10/J;
133: iJ11 = J00/J;
135: for (i = 0; i < NODES_PER_EL; i++) {
136: GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
137: GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
138: }
140: if (det_J) *det_J = J;
141: }
143: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
144: {
145: *ngp = 4;
146: gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919;
147: gp_xi[1][0] = -0.57735026919; gp_xi[1][1] = 0.57735026919;
148: gp_xi[2][0] = 0.57735026919; gp_xi[2][1] = 0.57735026919;
149: gp_xi[3][0] = 0.57735026919; gp_xi[3][1] = -0.57735026919;
150: gp_weight[0] = 1.0;
151: gp_weight[1] = 1.0;
152: gp_weight[2] = 1.0;
153: gp_weight[3] = 1.0;
154: }
156: /* procs to the left claim the ghost node as their element */
159: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
160: {
161: PetscInt m,n,p,M,N,P;
162: PetscInt sx,sy,sz;
165: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
166: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
168: if (mxl) {
169: *mxl = m;
170: if ((sx+m) == M) *mxl = m-1; /* last proc */
171: }
172: if (myl) {
173: *myl = n;
174: if ((sy+n) == N) *myl = n-1; /* last proc */
175: }
176: if (mzl) {
177: *mzl = p;
178: if ((sz+p) == P) *mzl = p-1; /* last proc */
179: }
180: return(0);
181: }
185: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
186: {
187: PetscInt si,sj,sk;
190: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
192: *sx = si;
193: if (si) *sx = si+1;
195: *sy = sj;
196: if (sj) *sy = sj+1;
198: if (sk) {
199: *sz = sk;
200: if (sk != 0) *sz = sk+1;
201: }
203: DMDAGetLocalElementSize(da,mx,my,mz);
204: return(0);
205: }
207: /*
208: i,j are the element indices
209: The unknown is a vector quantity.
210: The s[].c is used to indicate the degree of freedom.
211: */
214: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
215: {
217: /* velocity */
218: /* node 0 */
219: s_u[0].i = i; s_u[0].j = j; s_u[0].c = 0; /* Vx0 */
220: s_u[1].i = i; s_u[1].j = j; s_u[1].c = 1; /* Vy0 */
222: /* node 1 */
223: s_u[2].i = i; s_u[2].j = j+1; s_u[2].c = 0; /* Vx1 */
224: s_u[3].i = i; s_u[3].j = j+1; s_u[3].c = 1; /* Vy1 */
226: /* node 2 */
227: s_u[4].i = i+1; s_u[4].j = j+1; s_u[4].c = 0; /* Vx2 */
228: s_u[5].i = i+1; s_u[5].j = j+1; s_u[5].c = 1; /* Vy2 */
230: /* node 3 */
231: s_u[6].i = i+1; s_u[6].j = j; s_u[6].c = 0; /* Vx3 */
232: s_u[7].i = i+1; s_u[7].j = j; s_u[7].c = 1; /* Vy3 */
234: /* pressure */
235: s_p[0].i = i; s_p[0].j = j; s_p[0].c = 2; /* P0 */
236: s_p[1].i = i; s_p[1].j = j+1; s_p[1].c = 2; /* P1 */
237: s_p[2].i = i+1; s_p[2].j = j+1; s_p[2].c = 2; /* P2 */
238: s_p[3].i = i+1; s_p[3].j = j; s_p[3].c = 2; /* P3 */
239: return(0);
240: }
244: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
245: {
247: PetscMPIInt rank;
248: PetscInt proc_I,proc_J;
249: PetscInt cpu_x,cpu_y;
250: PetscInt local_mx,local_my;
251: Vec vlx,vly;
252: PetscInt *LX,*LY,i;
253: PetscScalar *_a;
254: Vec V_SEQ;
255: VecScatter ctx;
258: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
260: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
262: proc_J = rank/cpu_x;
263: proc_I = rank-cpu_x*proc_J;
265: PetscMalloc1(cpu_x,&LX);
266: PetscMalloc1(cpu_y,&LY);
268: DMDAGetLocalElementSize(da,&local_mx,&local_my,NULL);
269: VecCreate(PETSC_COMM_WORLD,&vlx);
270: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
271: VecSetFromOptions(vlx);
273: VecCreate(PETSC_COMM_WORLD,&vly);
274: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
275: VecSetFromOptions(vly);
277: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
278: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
279: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
280: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
282: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
283: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
284: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
285: VecGetArray(V_SEQ,&_a);
286: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
287: VecRestoreArray(V_SEQ,&_a);
288: VecScatterDestroy(&ctx);
289: VecDestroy(&V_SEQ);
291: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
292: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
293: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
294: VecGetArray(V_SEQ,&_a);
295: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
296: VecRestoreArray(V_SEQ,&_a);
297: VecScatterDestroy(&ctx);
298: VecDestroy(&V_SEQ);
300: *_lx = LX;
301: *_ly = LY;
303: VecDestroy(&vlx);
304: VecDestroy(&vly);
305: return(0);
306: }
310: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
311: {
312: DM cda;
313: Vec coords;
314: DMDACoor2d **_coords;
315: PetscInt si,sj,nx,ny,i,j;
316: FILE *fp;
317: char fname[PETSC_MAX_PATH_LEN];
318: PetscMPIInt rank;
322: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
323: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
324: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
325: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
327: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
329: DMGetCoordinateDM(da,&cda);
330: DMGetCoordinatesLocal(da,&coords);
331: DMDAVecGetArrayRead(cda,coords,&_coords);
332: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
333: for (j = sj; j < sj+ny-1; j++) {
334: for (i = si; i < si+nx-1; i++) {
335: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
336: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i].x),(double)PetscRealPart(_coords[j+1][i].y));
337: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i+1].x),(double)PetscRealPart(_coords[j+1][i+1].y));
338: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i+1].x),(double)PetscRealPart(_coords[j][i+1].y));
339: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
340: }
341: }
342: DMDAVecRestoreArrayRead(cda,coords,&_coords);
344: PetscFClose(PETSC_COMM_SELF,fp);
345: return(0);
346: }
350: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
351: {
352: DM cda;
353: Vec coords,local_fields;
354: DMDACoor2d **_coords;
355: FILE *fp;
356: char fname[PETSC_MAX_PATH_LEN];
357: PetscMPIInt rank;
358: PetscInt si,sj,nx,ny,i,j;
359: PetscInt n_dofs,d;
360: PetscScalar *_fields;
364: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
365: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
366: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
367: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
369: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
370: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
371: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
372: for (d = 0; d < n_dofs; d++) {
373: const char *field_name;
374: DMDAGetFieldName(da,d,&field_name);
375: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
376: }
377: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
379: DMGetCoordinateDM(da,&cda);
380: DMGetCoordinatesLocal(da,&coords);
381: DMDAVecGetArray(cda,coords,&_coords);
382: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
384: DMCreateLocalVector(da,&local_fields);
385: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
386: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
387: VecGetArray(local_fields,&_fields);
389: for (j = sj; j < sj+ny; j++) {
390: for (i = si; i < si+nx; i++) {
391: PetscScalar coord_x,coord_y;
392: PetscScalar field_d;
394: coord_x = _coords[j][i].x;
395: coord_y = _coords[j][i].y;
397: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
398: for (d = 0; d < n_dofs; d++) {
399: field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
400: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
401: }
402: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
403: }
404: }
405: VecRestoreArray(local_fields,&_fields);
406: VecDestroy(&local_fields);
408: DMDAVecRestoreArray(cda,coords,&_coords);
410: PetscFClose(PETSC_COMM_SELF,fp);
411: return(0);
412: }
416: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
417: {
418: DM cda;
419: Vec local_fields;
420: FILE *fp;
421: char fname[PETSC_MAX_PATH_LEN];
422: PetscMPIInt rank;
423: PetscInt si,sj,nx,ny,i,j,p;
424: PetscInt n_dofs,d;
425: GaussPointCoefficients **_coefficients;
426: PetscErrorCode ierr;
429: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
430: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
431: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
432: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
434: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
435: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
436: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
437: for (d = 0; d < n_dofs; d++) {
438: const char *field_name;
439: DMDAGetFieldName(da,d,&field_name);
440: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
441: }
442: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
444: DMGetCoordinateDM(da,&cda);
445: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
447: DMCreateLocalVector(da,&local_fields);
448: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
449: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
450: DMDAVecGetArray(da,local_fields,&_coefficients);
452: for (j = sj; j < sj+ny; j++) {
453: for (i = si; i < si+nx; i++) {
454: PetscScalar coord_x,coord_y;
456: for (p = 0; p < GAUSS_POINTS; p++) {
457: coord_x = _coefficients[j][i].gp_coords[2*p];
458: coord_y = _coefficients[j][i].gp_coords[2*p+1];
460: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
462: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e",(double)PetscRealPart(_coefficients[j][i].eta[p]),(double)PetscRealPart(_coefficients[j][i].fx[p]),(double)PetscRealPart(_coefficients[j][i].fy[p]));
463: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
464: }
465: }
466: }
467: DMDAVecRestoreArray(da,local_fields,&_coefficients);
468: VecDestroy(&local_fields);
470: PetscFClose(PETSC_COMM_SELF,fp);
471: return(0);
472: }
474: 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)
475: {
476: PetscInt ij;
477: PetscInt r,c,nc;
479: nc = u_NPE*u_dof;
481: r = w_dof*wi+wd;
482: c = u_dof*ui+ud;
484: ij = r*nc+c;
486: return ij;
487: }
489: /*
490: D = [ 2.eta 0 0 ]
491: [ 0 2.eta 0 ]
492: [ 0 0 eta ]
493:
494: B = [ d_dx 0 ]
495: [ 0 d_dy ]
496: [ d_dy d_dx ]
497: */
498: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
499: {
500: PetscInt ngp;
501: PetscScalar gp_xi[GAUSS_POINTS][2];
502: PetscScalar gp_weight[GAUSS_POINTS];
503: PetscInt p,i,j,k;
504: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
505: PetscScalar J_p,tildeD[3];
506: PetscScalar B[3][U_DOFS*NODES_PER_EL];
508: /* define quadrature rule */
509: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
511: /* evaluate integral */
512: for (p = 0; p < ngp; p++) {
513: ConstructQ12D_GNi(gp_xi[p],GNi_p);
514: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
516: for (i = 0; i < NODES_PER_EL; i++) {
517: PetscScalar d_dx_i = GNx_p[0][i];
518: PetscScalar d_dy_i = GNx_p[1][i];
520: B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
521: B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
522: B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
523: }
525: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
526: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
527: tildeD[2] = gp_weight[p]*J_p*eta[p];
529: /* form Bt tildeD B */
530: /*
531: Ke_ij = Bt_ik . D_kl . B_lj
532: = B_ki . D_kl . B_lj
533: = B_ki . D_kk . B_kj
534: */
535: for (i = 0; i < 8; i++) {
536: for (j = 0; j < 8; j++) {
537: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
538: Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
539: }
540: }
541: }
542: }
543: }
545: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
546: {
547: PetscInt ngp;
548: PetscScalar gp_xi[GAUSS_POINTS][2];
549: PetscScalar gp_weight[GAUSS_POINTS];
550: PetscInt p,i,j,di;
551: PetscScalar Ni_p[NODES_PER_EL];
552: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
553: PetscScalar J_p,fac;
555: /* define quadrature rule */
556: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
558: /* evaluate integral */
559: for (p = 0; p < ngp; p++) {
560: ConstructQ12D_Ni(gp_xi[p],Ni_p);
561: ConstructQ12D_GNi(gp_xi[p],GNi_p);
562: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
563: fac = gp_weight[p]*J_p;
565: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
566: for (di = 0; di < NSD; di++) { /* u dofs */
567: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
568: PetscInt IJ;
569: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);
571: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
572: }
573: }
574: }
575: }
576: }
578: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
579: {
580: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
581: PetscInt i,j;
582: PetscInt nr_g,nc_g;
584: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
585: FormGradientOperatorQ1(Ge,coords);
587: nr_g = U_DOFS*NODES_PER_EL;
588: nc_g = P_DOFS*NODES_PER_EL;
590: for (i = 0; i < nr_g; i++) {
591: for (j = 0; j < nc_g; j++) {
592: De[nr_g*j+i] = Ge[nc_g*i+j];
593: }
594: }
595: }
597: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
598: {
599: PetscInt ngp;
600: PetscScalar gp_xi[GAUSS_POINTS][2];
601: PetscScalar gp_weight[GAUSS_POINTS];
602: PetscInt p,i,j;
603: PetscScalar Ni_p[NODES_PER_EL];
604: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
605: PetscScalar J_p,fac,eta_avg;
607: /* define quadrature rule */
608: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
610: /* evaluate integral */
611: for (p = 0; p < ngp; p++) {
612: ConstructQ12D_Ni(gp_xi[p],Ni_p);
613: ConstructQ12D_GNi(gp_xi[p],GNi_p);
614: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
615: fac = gp_weight[p]*J_p;
617: for (i = 0; i < NODES_PER_EL; i++) {
618: for (j = 0; j < NODES_PER_EL; j++) {
619: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
620: }
621: }
622: }
624: /* scale */
625: eta_avg = 0.0;
626: for (p = 0; p < ngp; p++) eta_avg += eta[p];
627: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
628: fac = 1.0/eta_avg;
629: for (i = 0; i < NODES_PER_EL; i++) {
630: for (j = 0; j < NODES_PER_EL; j++) {
631: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
632: }
633: }
634: }
636: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
637: {
638: PetscInt ngp;
639: PetscScalar gp_xi[GAUSS_POINTS][2];
640: PetscScalar gp_weight[GAUSS_POINTS];
641: PetscInt p,i,j;
642: PetscScalar Ni_p[NODES_PER_EL];
643: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
644: PetscScalar J_p,fac,eta_avg;
646: /* define quadrature rule */
647: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
649: /* evaluate integral */
650: for (p = 0; p < ngp; p++) {
651: ConstructQ12D_Ni(gp_xi[p],Ni_p);
652: ConstructQ12D_GNi(gp_xi[p],GNi_p);
653: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
654: fac = gp_weight[p]*J_p;
656: for (i = 0; i < NODES_PER_EL; i++) {
657: for (j = 0; j < NODES_PER_EL; j++) {
658: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
659: }
660: }
661: }
663: /* scale */
664: eta_avg = 0.0;
665: for (p = 0; p < ngp; p++) eta_avg += eta[p];
666: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
667: fac = 1.0/eta_avg;
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] = fac*Ke[NODES_PER_EL*i+j];
671: }
672: }
673: }
675: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
676: {
677: PetscInt ngp;
678: PetscScalar gp_xi[GAUSS_POINTS][2];
679: PetscScalar gp_weight[GAUSS_POINTS];
680: PetscInt p,i;
681: PetscScalar Ni_p[NODES_PER_EL];
682: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
683: PetscScalar J_p,fac;
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: Fe[NSD*i] += fac*Ni_p[i]*fx[p];
697: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
698: }
699: }
700: }
704: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
705: {
707: /* get coords for the element */
708: el_coords[NSD*0] = _coords[ej][ei].x; el_coords[NSD*0+1] = _coords[ej][ei].y;
709: el_coords[NSD*1] = _coords[ej+1][ei].x; el_coords[NSD*1+1] = _coords[ej+1][ei].y;
710: el_coords[NSD*2] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
711: el_coords[NSD*3] = _coords[ej][ei+1].x; el_coords[NSD*3+1] = _coords[ej][ei+1].y;
712: return(0);
713: }
717: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
718: {
719: DM cda;
720: Vec coords;
721: DMDACoor2d **_coords;
722: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
723: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
724: PetscInt sex,sey,mx,my;
725: PetscInt ei,ej;
726: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
727: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
728: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
729: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
730: PetscScalar el_coords[NODES_PER_EL*NSD];
731: Vec local_properties;
732: GaussPointCoefficients **props;
733: PetscScalar *prop_eta;
734: PetscErrorCode ierr;
737: /* setup for coords */
738: DMGetCoordinateDM(stokes_da,&cda);
739: DMGetCoordinatesLocal(stokes_da,&coords);
740: DMDAVecGetArray(cda,coords,&_coords);
742: /* setup for coefficients */
743: DMCreateLocalVector(properties_da,&local_properties);
744: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
745: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
746: DMDAVecGetArray(properties_da,local_properties,&props);
748: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
749: for (ej = sey; ej < sey+my; ej++) {
750: for (ei = sex; ei < sex+mx; ei++) {
751: /* get coords for the element */
752: GetElementCoords(_coords,ei,ej,el_coords);
754: /* get coefficients for the element */
755: prop_eta = props[ej][ei].eta;
757: /* initialise element stiffness matrix */
758: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
759: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
760: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
761: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
763: /* form element stiffness matrix */
764: FormStressOperatorQ1(Ae,el_coords,prop_eta);
765: FormGradientOperatorQ1(Ge,el_coords);
766: FormDivergenceOperatorQ1(De,el_coords);
767: FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);
769: /* insert element matrix into global matrix */
770: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
771: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
772: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
773: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
774: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
775: }
776: }
777: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
778: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
780: DMDAVecRestoreArray(cda,coords,&_coords);
782: DMDAVecRestoreArray(properties_da,local_properties,&props);
783: VecDestroy(&local_properties);
784: return(0);
785: }
789: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
790: {
791: DM cda;
792: Vec coords;
793: DMDACoor2d **_coords;
794: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
795: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
796: PetscInt sex,sey,mx,my;
797: PetscInt ei,ej;
798: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
799: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
800: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
801: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
802: PetscScalar el_coords[NODES_PER_EL*NSD];
803: Vec local_properties;
804: GaussPointCoefficients **props;
805: PetscScalar *prop_eta;
806: PetscErrorCode ierr;
809: /* setup for coords */
810: DMGetCoordinateDM(stokes_da,&cda);
811: DMGetCoordinatesLocal(stokes_da,&coords);
812: DMDAVecGetArray(cda,coords,&_coords);
814: /* setup for coefficients */
815: DMCreateLocalVector(properties_da,&local_properties);
816: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
817: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
818: DMDAVecGetArray(properties_da,local_properties,&props);
820: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
821: for (ej = sey; ej < sey+my; ej++) {
822: for (ei = sex; ei < sex+mx; ei++) {
823: /* get coords for the element */
824: GetElementCoords(_coords,ei,ej,el_coords);
826: /* get coefficients for the element */
827: prop_eta = props[ej][ei].eta;
829: /* initialise element stiffness matrix */
830: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
831: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
832: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
833: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
835: /* form element stiffness matrix */
836: FormStressOperatorQ1(Ae,el_coords,prop_eta);
837: FormGradientOperatorQ1(Ge,el_coords);
838: FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);
840: /* insert element matrix into global matrix */
841: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
842: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
843: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
844: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
845: }
846: }
847: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
848: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
850: DMDAVecRestoreArray(cda,coords,&_coords);
852: DMDAVecRestoreArray(properties_da,local_properties,&props);
853: VecDestroy(&local_properties);
854: return(0);
855: }
859: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
860: {
861: PetscInt n;
864: for (n = 0; n < 4; n++) {
865: 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];
866: 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];
867: 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];
868: }
869: return(0);
870: }
874: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
875: {
876: DM cda;
877: Vec coords;
878: DMDACoor2d **_coords;
879: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
880: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
881: PetscInt sex,sey,mx,my;
882: PetscInt ei,ej;
883: PetscScalar Fe[NODES_PER_EL*U_DOFS];
884: PetscScalar He[NODES_PER_EL*P_DOFS];
885: PetscScalar el_coords[NODES_PER_EL*NSD];
886: Vec local_properties;
887: GaussPointCoefficients **props;
888: PetscScalar *prop_fx,*prop_fy;
889: Vec local_F;
890: StokesDOF **ff;
891: PetscErrorCode ierr;
894: /* setup for coords */
895: DMGetCoordinateDM(stokes_da,&cda);
896: DMGetCoordinatesLocal(stokes_da,&coords);
897: DMDAVecGetArray(cda,coords,&_coords);
899: /* setup for coefficients */
900: DMGetLocalVector(properties_da,&local_properties);
901: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
902: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
903: DMDAVecGetArray(properties_da,local_properties,&props);
905: /* get acces to the vector */
906: DMGetLocalVector(stokes_da,&local_F);
907: VecZeroEntries(local_F);
908: DMDAVecGetArray(stokes_da,local_F,&ff);
910: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
911: for (ej = sey; ej < sey+my; ej++) {
912: for (ei = sex; ei < sex+mx; ei++) {
913: /* get coords for the element */
914: GetElementCoords(_coords,ei,ej,el_coords);
916: /* get coefficients for the element */
917: prop_fx = props[ej][ei].fx;
918: prop_fy = props[ej][ei].fy;
920: /* initialise element stiffness matrix */
921: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
922: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
924: /* form element stiffness matrix */
925: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
927: /* insert element matrix into global matrix */
928: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
930: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
931: }
932: }
934: DMDAVecRestoreArray(stokes_da,local_F,&ff);
935: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
936: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
937: DMRestoreLocalVector(stokes_da,&local_F);
939: DMDAVecRestoreArray(cda,coords,&_coords);
941: DMDAVecRestoreArray(properties_da,local_properties,&props);
942: DMRestoreLocalVector(properties_da,&local_properties);
943: return(0);
944: }
948: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM *_da,Vec *_X)
949: {
950: DM da,cda;
951: Vec X;
952: StokesDOF **_stokes;
953: Vec coords;
954: DMDACoor2d **_coords;
955: PetscInt si,sj,ei,ej,i,j;
959: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
960: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,&da);
961: DMDASetFieldName(da,0,"anlytic_Vx");
962: DMDASetFieldName(da,1,"anlytic_Vy");
963: DMDASetFieldName(da,2,"analytic_P");
965: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.,0.);
967: DMGetCoordinatesLocal(da,&coords);
968: DMGetCoordinateDM(da,&cda);
969: DMDAVecGetArray(cda,coords,&_coords);
971: DMCreateGlobalVector(da,&X);
972: DMDAVecGetArray(da,X,&_stokes);
974: DMDAGetCorners(da,&si,&sj,0,&ei,&ej,0);
975: for (j = sj; j < sj+ej; j++) {
976: for (i = si; i < si+ei; i++) {
977: PetscReal pos[2],pressure,vel[2],total_stress[3],strain_rate[3];
979: pos[0] = PetscRealPart(_coords[j][i].x);
980: pos[1] = PetscRealPart(_coords[j][i].y);
982: evaluate_solCx(pos,eta0,eta1,xc,nz,vel,&pressure,total_stress,strain_rate);
984: _stokes[j][i].u_dof = vel[0];
985: _stokes[j][i].v_dof = vel[1];
986: _stokes[j][i].p_dof = pressure;
987: }
988: }
989: DMDAVecRestoreArray(da,X,&_stokes);
990: DMDAVecRestoreArray(cda,coords,&_coords);
992: *_da = da;
993: *_X = X;
994: return(0);
995: }
999: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
1000: {
1002: /* get the nodal fields */
1003: 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;
1004: 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;
1005: 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;
1006: 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;
1007: return(0);
1008: }
1012: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
1013: {
1014: DM cda;
1015: Vec coords,X_analytic_local,X_local;
1016: DMDACoor2d **_coords;
1017: PetscInt sex,sey,mx,my;
1018: PetscInt ei,ej;
1019: PetscScalar el_coords[NODES_PER_EL*NSD];
1020: StokesDOF **stokes_analytic,**stokes;
1021: StokesDOF stokes_analytic_e[4],stokes_e[4];
1023: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1024: PetscScalar Ni_p[NODES_PER_EL];
1025: PetscInt ngp;
1026: PetscScalar gp_xi[GAUSS_POINTS][2];
1027: PetscScalar gp_weight[GAUSS_POINTS];
1028: PetscInt p,i;
1029: PetscScalar J_p,fac;
1030: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1031: PetscInt M;
1032: PetscReal xymin[2],xymax[2];
1036: /* define quadrature rule */
1037: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1039: /* setup for coords */
1040: DMGetCoordinateDM(stokes_da,&cda);
1041: DMGetCoordinatesLocal(stokes_da,&coords);
1042: DMDAVecGetArray(cda,coords,&_coords);
1044: /* setup for analytic */
1045: DMCreateLocalVector(stokes_da,&X_analytic_local);
1046: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1047: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1048: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1050: /* setup for solution */
1051: DMCreateLocalVector(stokes_da,&X_local);
1052: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1053: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1054: DMDAVecGetArray(stokes_da,X_local,&stokes);
1056: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1057: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1059: h = (xymax[0]-xymin[0])/((PetscReal)M);
1061: tp_L2 = tu_L2 = tu_H1 = 0.0;
1063: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
1064: for (ej = sey; ej < sey+my; ej++) {
1065: for (ei = sex; ei < sex+mx; ei++) {
1066: /* get coords for the element */
1067: GetElementCoords(_coords,ei,ej,el_coords);
1068: StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1069: StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);
1071: /* evaluate integral */
1072: p_e_L2 = 0.0;
1073: u_e_L2 = 0.0;
1074: u_e_H1 = 0.0;
1075: for (p = 0; p < ngp; p++) {
1076: ConstructQ12D_Ni(gp_xi[p],Ni_p);
1077: ConstructQ12D_GNi(gp_xi[p],GNi_p);
1078: ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1079: fac = gp_weight[p]*J_p;
1081: for (i = 0; i < NODES_PER_EL; i++) {
1082: PetscScalar u_error,v_error;
1084: 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);
1086: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1087: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1088: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);
1090: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1091: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error /* du/dy */
1092: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1093: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error); /* dv/dy */
1094: }
1095: }
1097: tp_L2 += p_e_L2;
1098: tu_L2 += u_e_L2;
1099: tu_H1 += u_e_H1;
1100: }
1101: }
1102: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1103: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1104: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1105: p_L2 = PetscSqrtScalar(p_L2);
1106: u_L2 = PetscSqrtScalar(u_L2);
1107: u_H1 = PetscSqrtScalar(u_H1);
1109: PetscPrintf(PETSC_COMM_WORLD,"%1.4e %1.4e %1.4e %1.4e \n",(double)PetscRealPart(h),(double)PetscRealPart(p_L2),(double)PetscRealPart(u_L2),(double)PetscRealPart(u_H1));
1111: DMDAVecRestoreArray(cda,coords,&_coords);
1113: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1114: VecDestroy(&X_analytic_local);
1115: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1116: VecDestroy(&X_local);
1117: return(0);
1118: }
1122: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1123: {
1124: DM da_Stokes,da_prop;
1125: PetscInt u_dof,p_dof,dof,stencil_width;
1126: Mat A,B;
1127: PetscInt mxl,myl;
1128: DM prop_cda,vel_cda;
1129: Vec prop_coords,vel_coords;
1130: PetscInt si,sj,nx,ny,i,j,p;
1131: Vec f,X;
1132: PetscInt prop_dof,prop_stencil_width;
1133: Vec properties,l_properties;
1134: PetscReal dx,dy;
1135: PetscInt M,N;
1136: DMDACoor2d **_prop_coords,**_vel_coords;
1137: GaussPointCoefficients **element_props;
1138: PetscInt its;
1139: KSP ksp_S;
1140: PetscInt coefficient_structure = 0;
1141: PetscInt cpu_x,cpu_y,*lx = NULL,*ly = NULL;
1142: PetscBool use_gp_coords = PETSC_FALSE,set,output_gnuplot = PETSC_FALSE;
1143: char filename[PETSC_MAX_PATH_LEN];
1144: PetscErrorCode ierr;
1148: PetscOptionsGetBool(NULL,NULL,"-gnuplot",&output_gnuplot,NULL);
1149:
1150: /* Generate the da for velocity and pressure */
1151: /*
1152: We use Q1 elements for the temperature.
1153: FEM has a 9-point stencil (BOX) or connectivity pattern
1154: Num nodes in each direction is mx+1, my+1
1155: */
1156: u_dof = U_DOFS; /* Vx, Vy - velocities */
1157: p_dof = P_DOFS; /* p - pressure */
1158: dof = u_dof+p_dof;
1159: stencil_width = 1;
1160: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1161: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&da_Stokes);
1162: DMDASetFieldName(da_Stokes,0,"Vx");
1163: DMDASetFieldName(da_Stokes,1,"Vy");
1164: DMDASetFieldName(da_Stokes,2,"P");
1166: /* unit box [0,1] x [0,1] */
1167: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.,0.);
1169: /* Generate element properties, we will assume all material properties are constant over the element */
1170: /* local number of elements */
1171: DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,NULL);
1173: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1174: DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1175: DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);
1177: prop_dof = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1178: prop_stencil_width = 0;
1179: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1180: mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
1181: PetscFree(lx);
1182: PetscFree(ly);
1184: /* define centroid positions */
1185: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1186: dx = 1.0/((PetscReal)(M));
1187: dy = 1.0/((PetscReal)(N));
1189: 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);
1191: /* define coefficients */
1192: PetscOptionsGetInt(NULL,NULL,"-c_str",&coefficient_structure,NULL);
1194: DMCreateGlobalVector(da_prop,&properties);
1195: DMCreateLocalVector(da_prop,&l_properties);
1196: DMDAVecGetArray(da_prop,l_properties,&element_props);
1198: DMGetCoordinateDM(da_prop,&prop_cda);
1199: DMGetCoordinatesLocal(da_prop,&prop_coords);
1200: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
1202: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
1204: DMGetCoordinateDM(da_Stokes,&vel_cda);
1205: DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1206: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1208: /* interpolate the coordinates */
1209: for (j = sj; j < sj+ny; j++) {
1210: for (i = si; i < si+nx; i++) {
1211: PetscInt ngp;
1212: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1213: PetscScalar el_coords[8];
1215: GetElementCoords(_vel_coords,i,j,el_coords);
1216: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1218: for (p = 0; p < GAUSS_POINTS; p++) {
1219: PetscScalar gp_x,gp_y;
1220: PetscInt n;
1221: PetscScalar xi_p[2],Ni_p[4];
1223: xi_p[0] = gp_xi[p][0];
1224: xi_p[1] = gp_xi[p][1];
1225: ConstructQ12D_Ni(xi_p,Ni_p);
1227: gp_x = 0.0;
1228: gp_y = 0.0;
1229: for (n = 0; n < NODES_PER_EL; n++) {
1230: gp_x = gp_x+Ni_p[n]*el_coords[2*n];
1231: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1232: }
1233: element_props[j][i].gp_coords[2*p] = gp_x;
1234: element_props[j][i].gp_coords[2*p+1] = gp_y;
1235: }
1236: }
1237: }
1239: /* define the coefficients */
1240: PetscOptionsGetBool(NULL,NULL,"-use_gp_coords",&use_gp_coords,NULL);
1242: for (j = sj; j < sj+ny; j++) {
1243: for (i = si; i < si+nx; i++) {
1244: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1245: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1246: PetscReal coord_x,coord_y;
1248: if (coefficient_structure == 0) {
1249: PetscReal opts_eta0,opts_eta1,opts_xc;
1250: PetscInt opts_nz;
1252: opts_eta0 = 1.0;
1253: opts_eta1 = 1.0;
1254: opts_xc = 0.5;
1255: opts_nz = 1;
1257: PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1258: PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1259: PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1260: PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);
1262: for (p = 0; p < GAUSS_POINTS; p++) {
1263: coord_x = centroid_x;
1264: coord_y = centroid_y;
1265: if (use_gp_coords) {
1266: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1267: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1268: }
1270: element_props[j][i].eta[p] = opts_eta0;
1271: if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;
1273: element_props[j][i].fx[p] = 0.0;
1274: element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1275: }
1276: } else if (coefficient_structure == 1) { /* square sinker */
1277: PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;
1279: opts_eta0 = 1.0;
1280: opts_eta1 = 1.0;
1281: opts_dx = 0.50;
1282: opts_dy = 0.50;
1284: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1285: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1286: PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1287: PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
1289: for (p = 0; p < GAUSS_POINTS; p++) {
1290: coord_x = centroid_x;
1291: coord_y = centroid_y;
1292: if (use_gp_coords) {
1293: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1294: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1295: }
1297: element_props[j][i].eta[p] = opts_eta0;
1298: element_props[j][i].fx[p] = 0.0;
1299: element_props[j][i].fy[p] = 0.0;
1301: if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1302: if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1303: element_props[j][i].eta[p] = opts_eta1;
1304: element_props[j][i].fx[p] = 0.0;
1305: element_props[j][i].fy[p] = -1.0;
1306: }
1307: }
1308: }
1309: } else if (coefficient_structure == 2) { /* circular sinker */
1310: PetscReal opts_eta0,opts_eta1,opts_r,radius2;
1312: opts_eta0 = 1.0;
1313: opts_eta1 = 1.0;
1314: opts_r = 0.25;
1316: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1317: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1318: PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);
1320: for (p = 0; p < GAUSS_POINTS; p++) {
1321: coord_x = centroid_x;
1322: coord_y = centroid_y;
1323: if (use_gp_coords) {
1324: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1325: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1326: }
1328: element_props[j][i].eta[p] = opts_eta0;
1329: element_props[j][i].fx[p] = 0.0;
1330: element_props[j][i].fy[p] = 0.0;
1332: radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1333: if (radius2 < opts_r*opts_r) {
1334: element_props[j][i].eta[p] = opts_eta1;
1335: element_props[j][i].fx[p] = 0.0;
1336: element_props[j][i].fy[p] = -1.0;
1337: }
1338: }
1339: } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1340: PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;
1342: opts_eta0 = 1.0;
1343: opts_eta1 = 1.0;
1344: opts_r = 0.25;
1345: opts_c0x = 0.35; /* circle center */
1346: opts_c0y = 0.35;
1347: opts_s0x = 0.7; /* square center */
1348: opts_s0y = 0.7;
1349: opts_dx = 0.25;
1350: opts_dy = 0.25;
1351: opts_phi = 25;
1353: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1354: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1355: PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);
1356: PetscOptionsGetReal(NULL,NULL,"-sinker_c0x",&opts_c0x,NULL);
1357: PetscOptionsGetReal(NULL,NULL,"-sinker_c0y",&opts_c0y,NULL);
1358: PetscOptionsGetReal(NULL,NULL,"-sinker_s0x",&opts_s0x,NULL);
1359: PetscOptionsGetReal(NULL,NULL,"-sinker_s0y",&opts_s0y,NULL);
1360: PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1361: PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
1362: PetscOptionsGetReal(NULL,NULL,"-sinker_phi",&opts_phi,NULL);
1363: opts_phi *= PETSC_PI / 180;
1365: for (p = 0; p < GAUSS_POINTS; p++) {
1366: coord_x = centroid_x;
1367: coord_y = centroid_y;
1368: if (use_gp_coords) {
1369: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1370: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1371: }
1373: element_props[j][i].eta[p] = opts_eta0;
1374: element_props[j][i].fx[p] = 0.0;
1375: element_props[j][i].fy[p] = -0.2;
1377: radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1378: if (radius2 < opts_r*opts_r
1379: || (PetscAbs(+(coord_x - opts_s0x)*PetscCosReal(opts_phi) + (coord_y - opts_s0y)*PetscSinReal(opts_phi)) < opts_dx/2
1380: && PetscAbs(-(coord_x - opts_s0x)*PetscSinReal(opts_phi) + (coord_y - opts_s0y)*PetscCosReal(opts_phi)) < opts_dy/2)) {
1381: element_props[j][i].eta[p] = opts_eta1;
1382: element_props[j][i].fx[p] = 0.0;
1383: element_props[j][i].fy[p] = -1.0;
1384: }
1385: }
1386: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1387: }
1388: }
1389: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
1391: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1393: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1394: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1395: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
1397: if (output_gnuplot) {
1398: DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1399: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for Stokes eqn.","properties");
1400: }
1402: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1403: DMSetMatType(da_Stokes,MATAIJ);
1404: DMCreateMatrix(da_Stokes,&A);
1405: DMCreateMatrix(da_Stokes,&B);
1406: DMCreateGlobalVector(da_Stokes,&f);
1407: DMCreateGlobalVector(da_Stokes,&X);
1409: /* assemble A11 */
1410: MatZeroEntries(A);
1411: MatZeroEntries(B);
1412: VecZeroEntries(f);
1414: AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1415: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1416: /* build force vector */
1417: AssembleF_Stokes(f,da_Stokes,da_prop,properties);
1419: DMDABCApplyFreeSlip(da_Stokes,A,f);
1420: DMDABCApplyFreeSlip(da_Stokes,B,NULL);
1422: /* SOLVE */
1423: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1424: KSPSetOptionsPrefix(ksp_S,"stokes_");
1425: KSPSetDM(ksp_S,da_Stokes);
1426: KSPSetDMActive(ksp_S,PETSC_FALSE);
1427: KSPSetOperators(ksp_S,A,B);
1428: KSPSetFromOptions(ksp_S);
1429: {
1430: PC pc;
1431: const PetscInt ufields[] = {0,1},pfields[1] = {2};
1432: KSPGetPC(ksp_S,&pc);
1433: PCFieldSplitSetBlockSize(pc,3);
1434: PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1435: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1436: }
1438: {
1439: PC pc_S;
1440: KSP *sub_ksp,ksp_U;
1441: PetscInt nsplits;
1442: DM da_U;
1443: PetscBool is_pcfs;
1445: KSPSetUp(ksp_S);
1446: KSPGetPC(ksp_S,&pc_S);
1448: is_pcfs = PETSC_FALSE;
1449: PetscObjectTypeCompare((PetscObject)pc_S,PCFIELDSPLIT,&is_pcfs);
1451: if (is_pcfs) {
1452: PCFieldSplitGetSubKSP(pc_S,&nsplits,&sub_ksp);
1453: ksp_U = sub_ksp[0];
1454: PetscFree(sub_ksp);
1456: if (nsplits == 2) {
1457: DMDAGetReducedDMDA(da_Stokes,2,&da_U);
1459: KSPSetDM(ksp_U,da_U);
1460: KSPSetDMActive(ksp_U,PETSC_FALSE);
1461: DMDestroy(&da_U);
1462: }
1463: }
1464: }
1466: KSPSolve(ksp_S,f,X);
1468: PetscOptionsGetString(NULL,NULL,"-o",filename,sizeof(filename),&set);
1469: if (set) {
1470: char *ext;
1471: PetscViewer viewer;
1472: PetscBool flg;
1473: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1474: PetscStrrchr(filename,'.',&ext);
1475: PetscStrcmp("vts",ext,&flg);
1476: if (flg) {
1477: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1478: } else {
1479: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1480: }
1481: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1482: PetscViewerFileSetName(viewer,filename);
1483: VecView(X,viewer);
1484: PetscViewerDestroy(&viewer);
1485: }
1486: if (output_gnuplot) {
1487: DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
1488: }
1490: KSPGetIterationNumber(ksp_S,&its);
1492: if (coefficient_structure == 0) {
1493: PetscReal opts_eta0,opts_eta1,opts_xc;
1494: PetscInt opts_nz,N;
1495: DM da_Stokes_analytic;
1496: Vec X_analytic;
1497: PetscReal nrm1[3],nrm2[3],nrmI[3];
1499: opts_eta0 = 1.0;
1500: opts_eta1 = 1.0;
1501: opts_xc = 0.5;
1502: opts_nz = 1;
1504: PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1505: PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1506: PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1507: PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);
1509: DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1510: if (output_gnuplot) {
1511: DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
1512: }
1513: DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);
1515: VecAXPY(X_analytic,-1.0,X);
1516: VecGetSize(X_analytic,&N);
1517: N = N/3;
1519: VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1520: VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1521: VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);
1523: VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1524: VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1525: VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);
1527: VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1528: VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1529: VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);
1531: DMDestroy(&da_Stokes_analytic);
1532: VecDestroy(&X_analytic);
1533: }
1535: KSPDestroy(&ksp_S);
1536: VecDestroy(&X);
1537: VecDestroy(&f);
1538: MatDestroy(&A);
1539: MatDestroy(&B);
1541: DMDestroy(&da_Stokes);
1542: DMDestroy(&da_prop);
1544: VecDestroy(&properties);
1545: VecDestroy(&l_properties);
1546: return(0);
1547: }
1551: int main(int argc,char **args)
1552: {
1554: PetscInt mx,my;
1556: PetscInitialize(&argc,&args,(char*)0,help);
1558: mx = my = 10;
1559: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1560: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1562: solve_stokes_2d_coupled(mx,my);
1564: PetscFinalize();
1565: return 0;
1566: }
1568: /* -------------------------- helpers for boundary conditions -------------------------------- */
1571: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
1572: {
1573: DM cda;
1574: Vec coords;
1575: PetscInt si,sj,nx,ny,i,j;
1576: PetscInt M,N;
1577: DMDACoor2d **_coords;
1578: const PetscInt *g_idx;
1579: PetscInt *bc_global_ids;
1580: PetscScalar *bc_vals;
1581: PetscInt nbcs;
1582: PetscInt n_dofs;
1583: PetscErrorCode ierr;
1584: ISLocalToGlobalMapping ltogm;
1587: DMGetLocalToGlobalMapping(da,<ogm);
1588: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1590: DMGetCoordinateDM(da,&cda);
1591: DMGetCoordinatesLocal(da,&coords);
1592: DMDAVecGetArray(cda,coords,&_coords);
1593: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1594: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1596: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1597: PetscMalloc1(ny*n_dofs,&bc_vals);
1599: /* init the entries to -1 so VecSetValues will ignore them */
1600: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1602: i = nx-1;
1603: for (j = 0; j < ny; j++) {
1604: PetscInt local_id;
1606: local_id = i+j*nx;
1608: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1610: bc_vals[j] = 0.0;
1611: }
1612: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1613: nbcs = 0;
1614: if ((si+nx) == (M)) nbcs = ny;
1616: if (b) {
1617: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1618: VecAssemblyBegin(b);
1619: VecAssemblyEnd(b);
1620: }
1621: if (A) {
1622: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1623: }
1625: PetscFree(bc_vals);
1626: PetscFree(bc_global_ids);
1628: DMDAVecRestoreArray(cda,coords,&_coords);
1629: return(0);
1630: }
1634: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
1635: {
1636: DM cda;
1637: Vec coords;
1638: PetscInt si,sj,nx,ny,i,j;
1639: PetscInt M,N;
1640: DMDACoor2d **_coords;
1641: const PetscInt *g_idx;
1642: PetscInt *bc_global_ids;
1643: PetscScalar *bc_vals;
1644: PetscInt nbcs;
1645: PetscInt n_dofs;
1646: PetscErrorCode ierr;
1647: ISLocalToGlobalMapping ltogm;
1650: DMGetLocalToGlobalMapping(da,<ogm);
1651: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1653: DMGetCoordinateDM(da,&cda);
1654: DMGetCoordinatesLocal(da,&coords);
1655: DMDAVecGetArray(cda,coords,&_coords);
1656: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1657: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1659: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1660: PetscMalloc1(ny*n_dofs,&bc_vals);
1662: /* init the entries to -1 so VecSetValues will ignore them */
1663: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1665: i = 0;
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] = 0.0;
1674: }
1675: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1676: nbcs = 0;
1677: if (si == 0) nbcs = ny;
1679: if (b) {
1680: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1681: VecAssemblyBegin(b);
1682: VecAssemblyEnd(b);
1683: }
1685: if (A) {
1686: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1687: }
1689: PetscFree(bc_vals);
1690: PetscFree(bc_global_ids);
1692: DMDAVecRestoreArray(cda,coords,&_coords);
1693: return(0);
1694: }
1698: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
1699: {
1700: DM cda;
1701: Vec coords;
1702: PetscInt si,sj,nx,ny,i,j;
1703: PetscInt M,N;
1704: DMDACoor2d **_coords;
1705: const PetscInt *g_idx;
1706: PetscInt *bc_global_ids;
1707: PetscScalar *bc_vals;
1708: PetscInt nbcs;
1709: PetscInt n_dofs;
1710: PetscErrorCode ierr;
1711: ISLocalToGlobalMapping ltogm;
1714: DMGetLocalToGlobalMapping(da,<ogm);
1715: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1717: DMGetCoordinateDM(da,&cda);
1718: DMGetCoordinatesLocal(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: PetscMalloc1(nx,&bc_global_ids);
1724: PetscMalloc1(nx,&bc_vals);
1726: /* init the entries to -1 so VecSetValues will ignore them */
1727: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1729: j = ny-1;
1730: for (i = 0; i < nx; i++) {
1731: PetscInt local_id;
1733: local_id = i+j*nx;
1735: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1737: bc_vals[i] = 0.0;
1738: }
1739: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1740: nbcs = 0;
1741: if ((sj+ny) == (N)) nbcs = nx;
1743: if (b) {
1744: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1745: VecAssemblyBegin(b);
1746: VecAssemblyEnd(b);
1747: }
1748: if (A) {
1749: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1750: }
1752: PetscFree(bc_vals);
1753: PetscFree(bc_global_ids);
1755: DMDAVecRestoreArray(cda,coords,&_coords);
1756: return(0);
1757: }
1761: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
1762: {
1763: DM cda;
1764: Vec coords;
1765: PetscInt si,sj,nx,ny,i,j;
1766: PetscInt M,N;
1767: DMDACoor2d **_coords;
1768: const PetscInt *g_idx;
1769: PetscInt *bc_global_ids;
1770: PetscScalar *bc_vals;
1771: PetscInt nbcs;
1772: PetscInt n_dofs;
1773: PetscErrorCode ierr;
1774: ISLocalToGlobalMapping ltogm;
1777: DMGetLocalToGlobalMapping(da,<ogm);
1778: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1780: DMGetCoordinateDM(da,&cda);
1781: DMGetCoordinatesLocal(da,&coords);
1782: DMDAVecGetArray(cda,coords,&_coords);
1783: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1784: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1786: PetscMalloc1(nx,&bc_global_ids);
1787: PetscMalloc1(nx,&bc_vals);
1789: /* init the entries to -1 so VecSetValues will ignore them */
1790: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1792: j = 0;
1793: for (i = 0; i < nx; i++) {
1794: PetscInt local_id;
1796: local_id = i+j*nx;
1798: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1800: bc_vals[i] = 0.0;
1801: }
1802: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1803: nbcs = 0;
1804: if (sj == 0) nbcs = nx;
1806: if (b) {
1807: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1808: VecAssemblyBegin(b);
1809: VecAssemblyEnd(b);
1810: }
1811: if (A) {
1812: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1813: }
1815: PetscFree(bc_vals);
1816: PetscFree(bc_global_ids);
1818: DMDAVecRestoreArray(cda,coords,&_coords);
1819: return(0);
1820: }
1822: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
1825: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1826: {
1830: BCApplyZero_NORTH(da_Stokes,1,A,f);
1831: BCApplyZero_EAST(da_Stokes,0,A,f);
1832: BCApplyZero_SOUTH(da_Stokes,1,A,f);
1833: BCApplyZero_WEST(da_Stokes,0,A,f);
1834: return(0);
1835: }