Actual source code: ex43.c
petsc-3.8.4 2018-03-24
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 */
157: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
158: {
159: PetscInt m,n,p,M,N,P;
160: PetscInt sx,sy,sz;
163: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
164: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
166: if (mxl) {
167: *mxl = m;
168: if ((sx+m) == M) *mxl = m-1; /* last proc */
169: }
170: if (myl) {
171: *myl = n;
172: if ((sy+n) == N) *myl = n-1; /* last proc */
173: }
174: if (mzl) {
175: *mzl = p;
176: if ((sz+p) == P) *mzl = p-1; /* last proc */
177: }
178: return(0);
179: }
181: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
182: {
183: PetscInt si,sj,sk;
187: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
189: *sx = si;
190: if (si) *sx = si+1;
192: *sy = sj;
193: if (sj) *sy = sj+1;
195: if (sz) {
196: *sz = sk;
197: if (sk) *sz = sk+1;
198: }
200: DMDAGetLocalElementSize(da,mx,my,mz);
201: return(0);
202: }
204: /*
205: i,j are the element indices
206: The unknown is a vector quantity.
207: The s[].c is used to indicate the degree of freedom.
208: */
209: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
210: {
212: /* velocity */
213: /* node 0 */
214: s_u[0].i = i; s_u[0].j = j; s_u[0].c = 0; /* Vx0 */
215: s_u[1].i = i; s_u[1].j = j; s_u[1].c = 1; /* Vy0 */
217: /* node 1 */
218: s_u[2].i = i; s_u[2].j = j+1; s_u[2].c = 0; /* Vx1 */
219: s_u[3].i = i; s_u[3].j = j+1; s_u[3].c = 1; /* Vy1 */
221: /* node 2 */
222: s_u[4].i = i+1; s_u[4].j = j+1; s_u[4].c = 0; /* Vx2 */
223: s_u[5].i = i+1; s_u[5].j = j+1; s_u[5].c = 1; /* Vy2 */
225: /* node 3 */
226: s_u[6].i = i+1; s_u[6].j = j; s_u[6].c = 0; /* Vx3 */
227: s_u[7].i = i+1; s_u[7].j = j; s_u[7].c = 1; /* Vy3 */
229: /* pressure */
230: s_p[0].i = i; s_p[0].j = j; s_p[0].c = 2; /* P0 */
231: s_p[1].i = i; s_p[1].j = j+1; s_p[1].c = 2; /* P1 */
232: s_p[2].i = i+1; s_p[2].j = j+1; s_p[2].c = 2; /* P2 */
233: s_p[3].i = i+1; s_p[3].j = j; s_p[3].c = 2; /* P3 */
234: return(0);
235: }
237: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
238: {
240: PetscMPIInt rank;
241: PetscInt proc_I,proc_J;
242: PetscInt cpu_x,cpu_y;
243: PetscInt local_mx,local_my;
244: Vec vlx,vly;
245: PetscInt *LX,*LY,i;
246: PetscScalar *_a;
247: Vec V_SEQ;
248: VecScatter ctx;
251: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
253: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
255: proc_J = rank/cpu_x;
256: proc_I = rank-cpu_x*proc_J;
258: PetscMalloc1(cpu_x,&LX);
259: PetscMalloc1(cpu_y,&LY);
261: DMDAGetLocalElementSize(da,&local_mx,&local_my,NULL);
262: VecCreate(PETSC_COMM_WORLD,&vlx);
263: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
264: VecSetFromOptions(vlx);
266: VecCreate(PETSC_COMM_WORLD,&vly);
267: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
268: VecSetFromOptions(vly);
270: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
271: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
272: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
273: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
275: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
276: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
277: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
278: VecGetArray(V_SEQ,&_a);
279: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
280: VecRestoreArray(V_SEQ,&_a);
281: VecScatterDestroy(&ctx);
282: VecDestroy(&V_SEQ);
284: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
285: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
286: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
287: VecGetArray(V_SEQ,&_a);
288: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
289: VecRestoreArray(V_SEQ,&_a);
290: VecScatterDestroy(&ctx);
291: VecDestroy(&V_SEQ);
293: *_lx = LX;
294: *_ly = LY;
296: VecDestroy(&vlx);
297: VecDestroy(&vly);
298: return(0);
299: }
301: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
302: {
303: DM cda;
304: Vec coords;
305: DMDACoor2d **_coords;
306: PetscInt si,sj,nx,ny,i,j;
307: FILE *fp;
308: char fname[PETSC_MAX_PATH_LEN];
309: PetscMPIInt rank;
313: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
314: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
315: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
316: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
318: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
320: DMGetCoordinateDM(da,&cda);
321: DMGetCoordinatesLocal(da,&coords);
322: DMDAVecGetArrayRead(cda,coords,&_coords);
323: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
324: for (j = sj; j < sj+ny-1; j++) {
325: for (i = si; i < si+nx-1; i++) {
326: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
327: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i].x),(double)PetscRealPart(_coords[j+1][i].y));
328: 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));
329: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i+1].x),(double)PetscRealPart(_coords[j][i+1].y));
330: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
331: }
332: }
333: DMDAVecRestoreArrayRead(cda,coords,&_coords);
335: PetscFClose(PETSC_COMM_SELF,fp);
336: return(0);
337: }
339: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
340: {
341: DM cda;
342: Vec coords,local_fields;
343: DMDACoor2d **_coords;
344: FILE *fp;
345: char fname[PETSC_MAX_PATH_LEN];
346: PetscMPIInt rank;
347: PetscInt si,sj,nx,ny,i,j;
348: PetscInt n_dofs,d;
349: PetscScalar *_fields;
353: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
354: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
355: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
356: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
358: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
359: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
360: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
361: for (d = 0; d < n_dofs; d++) {
362: const char *field_name;
363: DMDAGetFieldName(da,d,&field_name);
364: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
365: }
366: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
368: DMGetCoordinateDM(da,&cda);
369: DMGetCoordinatesLocal(da,&coords);
370: DMDAVecGetArray(cda,coords,&_coords);
371: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
373: DMCreateLocalVector(da,&local_fields);
374: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
375: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
376: VecGetArray(local_fields,&_fields);
378: for (j = sj; j < sj+ny; j++) {
379: for (i = si; i < si+nx; i++) {
380: PetscScalar coord_x,coord_y;
381: PetscScalar field_d;
383: coord_x = _coords[j][i].x;
384: coord_y = _coords[j][i].y;
386: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
387: for (d = 0; d < n_dofs; d++) {
388: field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
389: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
390: }
391: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
392: }
393: }
394: VecRestoreArray(local_fields,&_fields);
395: VecDestroy(&local_fields);
397: DMDAVecRestoreArray(cda,coords,&_coords);
399: PetscFClose(PETSC_COMM_SELF,fp);
400: return(0);
401: }
403: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
404: {
405: DM cda;
406: Vec local_fields;
407: FILE *fp;
408: char fname[PETSC_MAX_PATH_LEN];
409: PetscMPIInt rank;
410: PetscInt si,sj,nx,ny,i,j,p;
411: PetscInt n_dofs,d;
412: GaussPointCoefficients **_coefficients;
413: PetscErrorCode ierr;
416: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
417: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
418: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
419: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
421: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
422: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
423: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
424: for (d = 0; d < n_dofs; d++) {
425: const char *field_name;
426: DMDAGetFieldName(da,d,&field_name);
427: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
428: }
429: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
431: DMGetCoordinateDM(da,&cda);
432: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
434: DMCreateLocalVector(da,&local_fields);
435: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
436: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
437: DMDAVecGetArray(da,local_fields,&_coefficients);
439: for (j = sj; j < sj+ny; j++) {
440: for (i = si; i < si+nx; i++) {
441: PetscScalar coord_x,coord_y;
443: for (p = 0; p < GAUSS_POINTS; p++) {
444: coord_x = _coefficients[j][i].gp_coords[2*p];
445: coord_y = _coefficients[j][i].gp_coords[2*p+1];
447: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
449: 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]));
450: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
451: }
452: }
453: }
454: DMDAVecRestoreArray(da,local_fields,&_coefficients);
455: VecDestroy(&local_fields);
457: PetscFClose(PETSC_COMM_SELF,fp);
458: return(0);
459: }
461: 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)
462: {
463: PetscInt ij;
464: PetscInt r,c,nc;
466: nc = u_NPE*u_dof;
467: r = w_dof*wi+wd;
468: c = u_dof*ui+ud;
469: ij = r*nc+c;
470: return ij;
471: }
473: /*
474: D = [ 2.eta 0 0 ]
475: [ 0 2.eta 0 ]
476: [ 0 0 eta ]
477:
478: B = [ d_dx 0 ]
479: [ 0 d_dy ]
480: [ d_dy d_dx ]
481: */
482: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
483: {
484: PetscInt ngp;
485: PetscScalar gp_xi[GAUSS_POINTS][2];
486: PetscScalar gp_weight[GAUSS_POINTS];
487: PetscInt p,i,j,k;
488: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
489: PetscScalar J_p,tildeD[3];
490: PetscScalar B[3][U_DOFS*NODES_PER_EL];
492: /* define quadrature rule */
493: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
495: /* evaluate integral */
496: for (p = 0; p < ngp; p++) {
497: ConstructQ12D_GNi(gp_xi[p],GNi_p);
498: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
500: for (i = 0; i < NODES_PER_EL; i++) {
501: PetscScalar d_dx_i = GNx_p[0][i];
502: PetscScalar d_dy_i = GNx_p[1][i];
504: B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
505: B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
506: B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
507: }
509: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
510: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
511: tildeD[2] = gp_weight[p]*J_p*eta[p];
513: /* form Bt tildeD B */
514: /*
515: Ke_ij = Bt_ik . D_kl . B_lj
516: = B_ki . D_kl . B_lj
517: = B_ki . D_kk . B_kj
518: */
519: for (i = 0; i < 8; i++) {
520: for (j = 0; j < 8; j++) {
521: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
522: Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
523: }
524: }
525: }
526: }
527: }
529: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
530: {
531: PetscInt ngp;
532: PetscScalar gp_xi[GAUSS_POINTS][2];
533: PetscScalar gp_weight[GAUSS_POINTS];
534: PetscInt p,i,j,di;
535: PetscScalar Ni_p[NODES_PER_EL];
536: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
537: PetscScalar J_p,fac;
539: /* define quadrature rule */
540: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
542: /* evaluate integral */
543: for (p = 0; p < ngp; p++) {
544: ConstructQ12D_Ni(gp_xi[p],Ni_p);
545: ConstructQ12D_GNi(gp_xi[p],GNi_p);
546: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
547: fac = gp_weight[p]*J_p;
549: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
550: for (di = 0; di < NSD; di++) { /* u dofs */
551: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
552: PetscInt IJ;
553: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);
555: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
556: }
557: }
558: }
559: }
560: }
562: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
563: {
564: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
565: PetscInt i,j;
566: PetscInt nr_g,nc_g;
568: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
569: FormGradientOperatorQ1(Ge,coords);
571: nr_g = U_DOFS*NODES_PER_EL;
572: nc_g = P_DOFS*NODES_PER_EL;
574: for (i = 0; i < nr_g; i++) {
575: for (j = 0; j < nc_g; j++) {
576: De[nr_g*j+i] = Ge[nc_g*i+j];
577: }
578: }
579: }
581: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
582: {
583: PetscInt ngp;
584: PetscScalar gp_xi[GAUSS_POINTS][2];
585: PetscScalar gp_weight[GAUSS_POINTS];
586: PetscInt p,i,j;
587: PetscScalar Ni_p[NODES_PER_EL];
588: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
589: PetscScalar J_p,fac,eta_avg;
591: /* define quadrature rule */
592: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
594: /* evaluate integral */
595: for (p = 0; p < ngp; p++) {
596: ConstructQ12D_Ni(gp_xi[p],Ni_p);
597: ConstructQ12D_GNi(gp_xi[p],GNi_p);
598: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
599: fac = gp_weight[p]*J_p;
601: for (i = 0; i < NODES_PER_EL; i++) {
602: for (j = 0; j < NODES_PER_EL; j++) {
603: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
604: }
605: }
606: }
608: /* scale */
609: eta_avg = 0.0;
610: for (p = 0; p < ngp; p++) eta_avg += eta[p];
611: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
612: fac = 1.0/eta_avg;
613: for (i = 0; i < NODES_PER_EL; i++) {
614: for (j = 0; j < NODES_PER_EL; j++) {
615: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
616: }
617: }
618: }
620: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
621: {
622: PetscInt ngp;
623: PetscScalar gp_xi[GAUSS_POINTS][2];
624: PetscScalar gp_weight[GAUSS_POINTS];
625: PetscInt p,i,j;
626: PetscScalar Ni_p[NODES_PER_EL];
627: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
628: PetscScalar J_p,fac,eta_avg;
630: /* define quadrature rule */
631: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
633: /* evaluate integral */
634: for (p = 0; p < ngp; p++) {
635: ConstructQ12D_Ni(gp_xi[p],Ni_p);
636: ConstructQ12D_GNi(gp_xi[p],GNi_p);
637: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
638: fac = gp_weight[p]*J_p;
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] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
643: }
644: }
645: }
647: /* scale */
648: eta_avg = 0.0;
649: for (p = 0; p < ngp; p++) eta_avg += eta[p];
650: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
651: fac = 1.0/eta_avg;
652: for (i = 0; i < NODES_PER_EL; i++) {
653: for (j = 0; j < NODES_PER_EL; j++) {
654: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
655: }
656: }
657: }
659: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
660: {
661: PetscInt ngp;
662: PetscScalar gp_xi[GAUSS_POINTS][2];
663: PetscScalar gp_weight[GAUSS_POINTS];
664: PetscInt p,i;
665: PetscScalar Ni_p[NODES_PER_EL];
666: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
667: PetscScalar J_p,fac;
669: /* define quadrature rule */
670: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
672: /* evaluate integral */
673: for (p = 0; p < ngp; p++) {
674: ConstructQ12D_Ni(gp_xi[p],Ni_p);
675: ConstructQ12D_GNi(gp_xi[p],GNi_p);
676: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
677: fac = gp_weight[p]*J_p;
679: for (i = 0; i < NODES_PER_EL; i++) {
680: Fe[NSD*i] += fac*Ni_p[i]*fx[p];
681: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
682: }
683: }
684: }
686: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
687: {
689: /* get coords for the element */
690: el_coords[NSD*0] = _coords[ej][ei].x; el_coords[NSD*0+1] = _coords[ej][ei].y;
691: el_coords[NSD*1] = _coords[ej+1][ei].x; el_coords[NSD*1+1] = _coords[ej+1][ei].y;
692: el_coords[NSD*2] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
693: el_coords[NSD*3] = _coords[ej][ei+1].x; el_coords[NSD*3+1] = _coords[ej][ei+1].y;
694: return(0);
695: }
697: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
698: {
699: DM cda;
700: Vec coords;
701: DMDACoor2d **_coords;
702: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
703: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
704: PetscInt sex,sey,mx,my;
705: PetscInt ei,ej;
706: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
707: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
708: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
709: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
710: PetscScalar el_coords[NODES_PER_EL*NSD];
711: Vec local_properties;
712: GaussPointCoefficients **props;
713: PetscScalar *prop_eta;
714: PetscErrorCode ierr;
717: /* setup for coords */
718: DMGetCoordinateDM(stokes_da,&cda);
719: DMGetCoordinatesLocal(stokes_da,&coords);
720: DMDAVecGetArray(cda,coords,&_coords);
722: /* setup for coefficients */
723: DMCreateLocalVector(properties_da,&local_properties);
724: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
725: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
726: DMDAVecGetArray(properties_da,local_properties,&props);
728: DMDAGetElementCorners(stokes_da,&sex,&sey,NULL,&mx,&my,NULL);
729: for (ej = sey; ej < sey+my; ej++) {
730: for (ei = sex; ei < sex+mx; ei++) {
731: /* get coords for the element */
732: GetElementCoords(_coords,ei,ej,el_coords);
734: /* get coefficients for the element */
735: prop_eta = props[ej][ei].eta;
737: /* initialise element stiffness matrix */
738: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
739: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
740: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
741: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
743: /* form element stiffness matrix */
744: FormStressOperatorQ1(Ae,el_coords,prop_eta);
745: FormGradientOperatorQ1(Ge,el_coords);
746: FormDivergenceOperatorQ1(De,el_coords);
747: FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);
749: /* insert element matrix into global matrix */
750: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
751: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
752: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
753: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
754: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
755: }
756: }
757: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
758: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
760: DMDAVecRestoreArray(cda,coords,&_coords);
762: DMDAVecRestoreArray(properties_da,local_properties,&props);
763: VecDestroy(&local_properties);
764: return(0);
765: }
767: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
768: {
769: DM cda;
770: Vec coords;
771: DMDACoor2d **_coords;
772: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
773: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
774: PetscInt sex,sey,mx,my;
775: PetscInt ei,ej;
776: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
777: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
778: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
779: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
780: PetscScalar el_coords[NODES_PER_EL*NSD];
781: Vec local_properties;
782: GaussPointCoefficients **props;
783: PetscScalar *prop_eta;
784: PetscErrorCode ierr;
787: /* setup for coords */
788: DMGetCoordinateDM(stokes_da,&cda);
789: DMGetCoordinatesLocal(stokes_da,&coords);
790: DMDAVecGetArray(cda,coords,&_coords);
792: /* setup for coefficients */
793: DMCreateLocalVector(properties_da,&local_properties);
794: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
795: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
796: DMDAVecGetArray(properties_da,local_properties,&props);
798: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
799: for (ej = sey; ej < sey+my; ej++) {
800: for (ei = sex; ei < sex+mx; ei++) {
801: /* get coords for the element */
802: GetElementCoords(_coords,ei,ej,el_coords);
804: /* get coefficients for the element */
805: prop_eta = props[ej][ei].eta;
807: /* initialise element stiffness matrix */
808: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
809: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
810: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
811: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
813: /* form element stiffness matrix */
814: FormStressOperatorQ1(Ae,el_coords,prop_eta);
815: FormGradientOperatorQ1(Ge,el_coords);
816: FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);
818: /* insert element matrix into global matrix */
819: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
820: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
821: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
822: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
823: }
824: }
825: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
826: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
828: DMDAVecRestoreArray(cda,coords,&_coords);
830: DMDAVecRestoreArray(properties_da,local_properties,&props);
831: VecDestroy(&local_properties);
832: return(0);
833: }
835: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
836: {
837: PetscInt n;
840: for (n = 0; n < 4; n++) {
841: 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];
842: 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];
843: 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];
844: }
845: return(0);
846: }
848: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
849: {
850: DM cda;
851: Vec coords;
852: DMDACoor2d **_coords;
853: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
854: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
855: PetscInt sex,sey,mx,my;
856: PetscInt ei,ej;
857: PetscScalar Fe[NODES_PER_EL*U_DOFS];
858: PetscScalar He[NODES_PER_EL*P_DOFS];
859: PetscScalar el_coords[NODES_PER_EL*NSD];
860: Vec local_properties;
861: GaussPointCoefficients **props;
862: PetscScalar *prop_fx,*prop_fy;
863: Vec local_F;
864: StokesDOF **ff;
865: PetscErrorCode ierr;
868: /* setup for coords */
869: DMGetCoordinateDM(stokes_da,&cda);
870: DMGetCoordinatesLocal(stokes_da,&coords);
871: DMDAVecGetArray(cda,coords,&_coords);
873: /* setup for coefficients */
874: DMGetLocalVector(properties_da,&local_properties);
875: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
876: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
877: DMDAVecGetArray(properties_da,local_properties,&props);
879: /* get acces to the vector */
880: DMGetLocalVector(stokes_da,&local_F);
881: VecZeroEntries(local_F);
882: DMDAVecGetArray(stokes_da,local_F,&ff);
884: DMDAGetElementCorners(stokes_da,&sex,&sey,NULL,&mx,&my,NULL);
885: for (ej = sey; ej < sey+my; ej++) {
886: for (ei = sex; ei < sex+mx; ei++) {
887: /* get coords for the element */
888: GetElementCoords(_coords,ei,ej,el_coords);
890: /* get coefficients for the element */
891: prop_fx = props[ej][ei].fx;
892: prop_fy = props[ej][ei].fy;
894: /* initialise element stiffness matrix */
895: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
896: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
898: /* form element stiffness matrix */
899: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
901: /* insert element matrix into global matrix */
902: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
904: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
905: }
906: }
908: DMDAVecRestoreArray(stokes_da,local_F,&ff);
909: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
910: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
911: DMRestoreLocalVector(stokes_da,&local_F);
913: DMDAVecRestoreArray(cda,coords,&_coords);
915: DMDAVecRestoreArray(properties_da,local_properties,&props);
916: DMRestoreLocalVector(properties_da,&local_properties);
917: return(0);
918: }
920: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM *_da,Vec *_X)
921: {
922: DM da,cda;
923: Vec X;
924: StokesDOF **_stokes;
925: Vec coords;
926: DMDACoor2d **_coords;
927: PetscInt si,sj,ei,ej,i,j;
931: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,3,1,NULL,NULL,&da);
932: DMSetFromOptions(da);
933: DMSetUp(da);
934: DMDASetFieldName(da,0,"anlytic_Vx");
935: DMDASetFieldName(da,1,"anlytic_Vy");
936: DMDASetFieldName(da,2,"analytic_P");
938: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.,0.);
940: DMGetCoordinatesLocal(da,&coords);
941: DMGetCoordinateDM(da,&cda);
942: DMDAVecGetArray(cda,coords,&_coords);
944: DMCreateGlobalVector(da,&X);
945: DMDAVecGetArray(da,X,&_stokes);
947: DMDAGetCorners(da,&si,&sj,0,&ei,&ej,0);
948: for (j = sj; j < sj+ej; j++) {
949: for (i = si; i < si+ei; i++) {
950: PetscReal pos[2],pressure,vel[2],total_stress[3],strain_rate[3];
952: pos[0] = PetscRealPart(_coords[j][i].x);
953: pos[1] = PetscRealPart(_coords[j][i].y);
955: evaluate_solCx(pos,eta0,eta1,xc,nz,vel,&pressure,total_stress,strain_rate);
957: _stokes[j][i].u_dof = vel[0];
958: _stokes[j][i].v_dof = vel[1];
959: _stokes[j][i].p_dof = pressure;
960: }
961: }
962: DMDAVecRestoreArray(da,X,&_stokes);
963: DMDAVecRestoreArray(cda,coords,&_coords);
965: *_da = da;
966: *_X = X;
967: return(0);
968: }
970: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
971: {
973: /* get the nodal fields */
974: 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;
975: 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;
976: 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;
977: 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;
978: return(0);
979: }
981: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
982: {
983: DM cda;
984: Vec coords,X_analytic_local,X_local;
985: DMDACoor2d **_coords;
986: PetscInt sex,sey,mx,my;
987: PetscInt ei,ej;
988: PetscScalar el_coords[NODES_PER_EL*NSD];
989: StokesDOF **stokes_analytic,**stokes;
990: StokesDOF stokes_analytic_e[4],stokes_e[4];
992: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
993: PetscScalar Ni_p[NODES_PER_EL];
994: PetscInt ngp;
995: PetscScalar gp_xi[GAUSS_POINTS][2];
996: PetscScalar gp_weight[GAUSS_POINTS];
997: PetscInt p,i;
998: PetscScalar J_p,fac;
999: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1000: PetscInt M;
1001: PetscReal xymin[2],xymax[2];
1005: /* define quadrature rule */
1006: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1008: /* setup for coords */
1009: DMGetCoordinateDM(stokes_da,&cda);
1010: DMGetCoordinatesLocal(stokes_da,&coords);
1011: DMDAVecGetArray(cda,coords,&_coords);
1013: /* setup for analytic */
1014: DMCreateLocalVector(stokes_da,&X_analytic_local);
1015: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1016: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1017: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1019: /* setup for solution */
1020: DMCreateLocalVector(stokes_da,&X_local);
1021: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1022: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1023: DMDAVecGetArray(stokes_da,X_local,&stokes);
1025: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1026: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1028: h = (xymax[0]-xymin[0])/((PetscReal)M);
1030: tp_L2 = tu_L2 = tu_H1 = 0.0;
1032: DMDAGetElementCorners(stokes_da,&sex,&sey,NULL,&mx,&my,NULL);
1033: for (ej = sey; ej < sey+my; ej++) {
1034: for (ei = sex; ei < sex+mx; ei++) {
1035: /* get coords for the element */
1036: GetElementCoords(_coords,ei,ej,el_coords);
1037: StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1038: StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);
1040: /* evaluate integral */
1041: p_e_L2 = 0.0;
1042: u_e_L2 = 0.0;
1043: u_e_H1 = 0.0;
1044: for (p = 0; p < ngp; p++) {
1045: ConstructQ12D_Ni(gp_xi[p],Ni_p);
1046: ConstructQ12D_GNi(gp_xi[p],GNi_p);
1047: ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1048: fac = gp_weight[p]*J_p;
1050: for (i = 0; i < NODES_PER_EL; i++) {
1051: PetscScalar u_error,v_error;
1053: 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);
1055: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1056: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1057: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);
1059: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1060: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error /* du/dy */
1061: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1062: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error); /* dv/dy */
1063: }
1064: }
1066: tp_L2 += p_e_L2;
1067: tu_L2 += u_e_L2;
1068: tu_H1 += u_e_H1;
1069: }
1070: }
1071: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1072: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1073: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1074: p_L2 = PetscSqrtScalar(p_L2);
1075: u_L2 = PetscSqrtScalar(u_L2);
1076: u_H1 = PetscSqrtScalar(u_H1);
1078: 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));
1080: DMDAVecRestoreArray(cda,coords,&_coords);
1082: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1083: VecDestroy(&X_analytic_local);
1084: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1085: VecDestroy(&X_local);
1086: return(0);
1087: }
1089: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1090: {
1091: DM da_Stokes,da_prop;
1092: PetscInt u_dof,p_dof,dof,stencil_width;
1093: Mat A,B;
1094: PetscInt mxl,myl;
1095: DM prop_cda,vel_cda;
1096: Vec prop_coords,vel_coords;
1097: PetscInt si,sj,nx,ny,i,j,p;
1098: Vec f,X;
1099: PetscInt prop_dof,prop_stencil_width;
1100: Vec properties,l_properties;
1101: PetscReal dx,dy;
1102: PetscInt M,N;
1103: DMDACoor2d **_prop_coords,**_vel_coords;
1104: GaussPointCoefficients **element_props;
1105: PetscInt its;
1106: KSP ksp_S;
1107: PetscInt coefficient_structure = 0;
1108: PetscInt cpu_x,cpu_y,*lx = NULL,*ly = NULL;
1109: PetscBool use_gp_coords = PETSC_FALSE,set,output_gnuplot = PETSC_FALSE;
1110: char filename[PETSC_MAX_PATH_LEN];
1111: PetscErrorCode ierr;
1115: PetscOptionsGetBool(NULL,NULL,"-gnuplot",&output_gnuplot,NULL);
1116:
1117: /* Generate the da for velocity and pressure */
1118: /*
1119: We use Q1 elements for the temperature.
1120: FEM has a 9-point stencil (BOX) or connectivity pattern
1121: Num nodes in each direction is mx+1, my+1
1122: */
1123: u_dof = U_DOFS; /* Vx, Vy - velocities */
1124: p_dof = P_DOFS; /* p - pressure */
1125: dof = u_dof+p_dof;
1126: stencil_width = 1;
1127: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&da_Stokes);
1128: DMSetFromOptions(da_Stokes);
1129: DMSetUp(da_Stokes);
1130: DMDASetFieldName(da_Stokes,0,"Vx");
1131: DMDASetFieldName(da_Stokes,1,"Vy");
1132: DMDASetFieldName(da_Stokes,2,"P");
1134: /* unit box [0,1] x [0,1] */
1135: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.,0.);
1137: /* Generate element properties, we will assume all material properties are constant over the element */
1138: /* local number of elements */
1139: DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,NULL);
1141: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1142: DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1143: DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);
1145: prop_dof = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1146: prop_stencil_width = 0;
1147: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
1148: DMSetFromOptions(da_prop);
1149: DMSetUp(da_prop);
1150: PetscFree(lx);
1151: PetscFree(ly);
1153: /* define centroid positions */
1154: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1155: dx = 1.0/((PetscReal)(M));
1156: dy = 1.0/((PetscReal)(N));
1158: 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);
1160: /* define coefficients */
1161: PetscOptionsGetInt(NULL,NULL,"-c_str",&coefficient_structure,NULL);
1163: DMCreateGlobalVector(da_prop,&properties);
1164: DMCreateLocalVector(da_prop,&l_properties);
1165: DMDAVecGetArray(da_prop,l_properties,&element_props);
1167: DMGetCoordinateDM(da_prop,&prop_cda);
1168: DMGetCoordinatesLocal(da_prop,&prop_coords);
1169: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
1171: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
1173: DMGetCoordinateDM(da_Stokes,&vel_cda);
1174: DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1175: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1177: /* interpolate the coordinates */
1178: for (j = sj; j < sj+ny; j++) {
1179: for (i = si; i < si+nx; i++) {
1180: PetscInt ngp;
1181: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1182: PetscScalar el_coords[8];
1184: GetElementCoords(_vel_coords,i,j,el_coords);
1185: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1187: for (p = 0; p < GAUSS_POINTS; p++) {
1188: PetscScalar gp_x,gp_y;
1189: PetscInt n;
1190: PetscScalar xi_p[2],Ni_p[4];
1192: xi_p[0] = gp_xi[p][0];
1193: xi_p[1] = gp_xi[p][1];
1194: ConstructQ12D_Ni(xi_p,Ni_p);
1196: gp_x = 0.0;
1197: gp_y = 0.0;
1198: for (n = 0; n < NODES_PER_EL; n++) {
1199: gp_x = gp_x+Ni_p[n]*el_coords[2*n];
1200: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1201: }
1202: element_props[j][i].gp_coords[2*p] = gp_x;
1203: element_props[j][i].gp_coords[2*p+1] = gp_y;
1204: }
1205: }
1206: }
1208: /* define the coefficients */
1209: PetscOptionsGetBool(NULL,NULL,"-use_gp_coords",&use_gp_coords,NULL);
1211: for (j = sj; j < sj+ny; j++) {
1212: for (i = si; i < si+nx; i++) {
1213: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1214: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1215: PetscReal coord_x,coord_y;
1217: if (coefficient_structure == 0) {
1218: PetscReal opts_eta0,opts_eta1,opts_xc;
1219: PetscInt opts_nz;
1221: opts_eta0 = 1.0;
1222: opts_eta1 = 1.0;
1223: opts_xc = 0.5;
1224: opts_nz = 1;
1226: PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1227: PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1228: PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1229: PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);
1231: for (p = 0; p < GAUSS_POINTS; p++) {
1232: coord_x = centroid_x;
1233: coord_y = centroid_y;
1234: if (use_gp_coords) {
1235: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1236: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1237: }
1239: element_props[j][i].eta[p] = opts_eta0;
1240: if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;
1242: element_props[j][i].fx[p] = 0.0;
1243: element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1244: }
1245: } else if (coefficient_structure == 1) { /* square sinker */
1246: PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;
1248: opts_eta0 = 1.0;
1249: opts_eta1 = 1.0;
1250: opts_dx = 0.50;
1251: opts_dy = 0.50;
1253: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1254: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1255: PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1256: PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
1258: for (p = 0; p < GAUSS_POINTS; p++) {
1259: coord_x = centroid_x;
1260: coord_y = centroid_y;
1261: if (use_gp_coords) {
1262: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1263: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1264: }
1266: element_props[j][i].eta[p] = opts_eta0;
1267: element_props[j][i].fx[p] = 0.0;
1268: element_props[j][i].fy[p] = 0.0;
1270: if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1271: if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1272: element_props[j][i].eta[p] = opts_eta1;
1273: element_props[j][i].fx[p] = 0.0;
1274: element_props[j][i].fy[p] = -1.0;
1275: }
1276: }
1277: }
1278: } else if (coefficient_structure == 2) { /* circular sinker */
1279: PetscReal opts_eta0,opts_eta1,opts_r,radius2;
1281: opts_eta0 = 1.0;
1282: opts_eta1 = 1.0;
1283: opts_r = 0.25;
1285: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1286: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1287: PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,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: radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1302: if (radius2 < opts_r*opts_r) {
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: } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1309: PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;
1311: opts_eta0 = 1.0;
1312: opts_eta1 = 1.0;
1313: opts_r = 0.25;
1314: opts_c0x = 0.35; /* circle center */
1315: opts_c0y = 0.35;
1316: opts_s0x = 0.7; /* square center */
1317: opts_s0y = 0.7;
1318: opts_dx = 0.25;
1319: opts_dy = 0.25;
1320: opts_phi = 25;
1322: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1323: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1324: PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);
1325: PetscOptionsGetReal(NULL,NULL,"-sinker_c0x",&opts_c0x,NULL);
1326: PetscOptionsGetReal(NULL,NULL,"-sinker_c0y",&opts_c0y,NULL);
1327: PetscOptionsGetReal(NULL,NULL,"-sinker_s0x",&opts_s0x,NULL);
1328: PetscOptionsGetReal(NULL,NULL,"-sinker_s0y",&opts_s0y,NULL);
1329: PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1330: PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
1331: PetscOptionsGetReal(NULL,NULL,"-sinker_phi",&opts_phi,NULL);
1332: opts_phi *= PETSC_PI / 180;
1334: for (p = 0; p < GAUSS_POINTS; p++) {
1335: coord_x = centroid_x;
1336: coord_y = centroid_y;
1337: if (use_gp_coords) {
1338: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1339: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1340: }
1342: element_props[j][i].eta[p] = opts_eta0;
1343: element_props[j][i].fx[p] = 0.0;
1344: element_props[j][i].fy[p] = -0.2;
1346: radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1347: if (radius2 < opts_r*opts_r
1348: || (PetscAbs(+(coord_x - opts_s0x)*PetscCosReal(opts_phi) + (coord_y - opts_s0y)*PetscSinReal(opts_phi)) < opts_dx/2
1349: && PetscAbs(-(coord_x - opts_s0x)*PetscSinReal(opts_phi) + (coord_y - opts_s0y)*PetscCosReal(opts_phi)) < opts_dy/2)) {
1350: element_props[j][i].eta[p] = opts_eta1;
1351: element_props[j][i].fx[p] = 0.0;
1352: element_props[j][i].fy[p] = -1.0;
1353: }
1354: }
1355: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1356: }
1357: }
1358: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
1360: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1362: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1363: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1364: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
1366: if (output_gnuplot) {
1367: DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1368: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for Stokes eqn.","properties");
1369: }
1371: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1372: DMSetMatType(da_Stokes,MATAIJ);
1373: DMCreateMatrix(da_Stokes,&A);
1374: DMCreateMatrix(da_Stokes,&B);
1375: DMCreateGlobalVector(da_Stokes,&f);
1376: DMCreateGlobalVector(da_Stokes,&X);
1378: /* assemble A11 */
1379: MatZeroEntries(A);
1380: MatZeroEntries(B);
1381: VecZeroEntries(f);
1383: AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1384: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1385: /* build force vector */
1386: AssembleF_Stokes(f,da_Stokes,da_prop,properties);
1388: DMDABCApplyFreeSlip(da_Stokes,A,f);
1389: DMDABCApplyFreeSlip(da_Stokes,B,NULL);
1391: /* SOLVE */
1392: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1393: KSPSetOptionsPrefix(ksp_S,"stokes_");
1394: KSPSetDM(ksp_S,da_Stokes);
1395: KSPSetDMActive(ksp_S,PETSC_FALSE);
1396: KSPSetOperators(ksp_S,A,B);
1397: KSPSetFromOptions(ksp_S);
1398: {
1399: PC pc;
1400: const PetscInt ufields[] = {0,1},pfields[1] = {2};
1401: KSPGetPC(ksp_S,&pc);
1402: PCFieldSplitSetBlockSize(pc,3);
1403: PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1404: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1405: }
1407: {
1408: PC pc_S;
1409: KSP *sub_ksp,ksp_U;
1410: PetscInt nsplits;
1411: DM da_U;
1412: PetscBool is_pcfs;
1414: KSPSetUp(ksp_S);
1415: KSPGetPC(ksp_S,&pc_S);
1417: is_pcfs = PETSC_FALSE;
1418: PetscObjectTypeCompare((PetscObject)pc_S,PCFIELDSPLIT,&is_pcfs);
1420: if (is_pcfs) {
1421: PCFieldSplitGetSubKSP(pc_S,&nsplits,&sub_ksp);
1422: ksp_U = sub_ksp[0];
1423: PetscFree(sub_ksp);
1425: if (nsplits == 2) {
1426: DMDAGetReducedDMDA(da_Stokes,2,&da_U);
1428: KSPSetDM(ksp_U,da_U);
1429: KSPSetDMActive(ksp_U,PETSC_FALSE);
1430: DMDestroy(&da_U);
1431: }
1432: }
1433: }
1435: KSPSolve(ksp_S,f,X);
1437: PetscOptionsGetString(NULL,NULL,"-o",filename,sizeof(filename),&set);
1438: if (set) {
1439: char *ext;
1440: PetscViewer viewer;
1441: PetscBool flg;
1442: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1443: PetscStrrchr(filename,'.',&ext);
1444: PetscStrcmp("vts",ext,&flg);
1445: if (flg) {
1446: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1447: } else {
1448: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1449: }
1450: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1451: PetscViewerFileSetName(viewer,filename);
1452: VecView(X,viewer);
1453: PetscViewerDestroy(&viewer);
1454: }
1455: if (output_gnuplot) {
1456: DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
1457: }
1459: KSPGetIterationNumber(ksp_S,&its);
1461: if (coefficient_structure == 0) {
1462: PetscReal opts_eta0,opts_eta1,opts_xc;
1463: PetscInt opts_nz,N;
1464: DM da_Stokes_analytic;
1465: Vec X_analytic;
1466: PetscReal nrm1[3],nrm2[3],nrmI[3];
1468: opts_eta0 = 1.0;
1469: opts_eta1 = 1.0;
1470: opts_xc = 0.5;
1471: opts_nz = 1;
1473: PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1474: PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1475: PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1476: PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);
1478: DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1479: if (output_gnuplot) {
1480: DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
1481: }
1482: DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);
1484: VecAXPY(X_analytic,-1.0,X);
1485: VecGetSize(X_analytic,&N);
1486: N = N/3;
1488: VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1489: VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1490: VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);
1492: VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1493: VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1494: VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);
1496: VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1497: VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1498: VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);
1500: DMDestroy(&da_Stokes_analytic);
1501: VecDestroy(&X_analytic);
1502: }
1504: KSPDestroy(&ksp_S);
1505: VecDestroy(&X);
1506: VecDestroy(&f);
1507: MatDestroy(&A);
1508: MatDestroy(&B);
1510: DMDestroy(&da_Stokes);
1511: DMDestroy(&da_prop);
1513: VecDestroy(&properties);
1514: VecDestroy(&l_properties);
1515: return(0);
1516: }
1518: int main(int argc,char **args)
1519: {
1521: PetscInt mx,my;
1523: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
1524: mx = my = 10;
1525: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1526: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1527: solve_stokes_2d_coupled(mx,my);
1528: PetscFinalize();
1529: return ierr;
1530: }
1532: /* -------------------------- helpers for boundary conditions -------------------------------- */
1533: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
1534: {
1535: DM cda;
1536: Vec coords;
1537: PetscInt si,sj,nx,ny,i,j;
1538: PetscInt M,N;
1539: DMDACoor2d **_coords;
1540: const PetscInt *g_idx;
1541: PetscInt *bc_global_ids;
1542: PetscScalar *bc_vals;
1543: PetscInt nbcs;
1544: PetscInt n_dofs;
1545: PetscErrorCode ierr;
1546: ISLocalToGlobalMapping ltogm;
1549: DMGetLocalToGlobalMapping(da,<ogm);
1550: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1552: DMGetCoordinateDM(da,&cda);
1553: DMGetCoordinatesLocal(da,&coords);
1554: DMDAVecGetArray(cda,coords,&_coords);
1555: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1556: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1558: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1559: PetscMalloc1(ny*n_dofs,&bc_vals);
1561: /* init the entries to -1 so VecSetValues will ignore them */
1562: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1564: i = nx-1;
1565: for (j = 0; j < ny; j++) {
1566: PetscInt local_id;
1568: local_id = i+j*nx;
1570: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1572: bc_vals[j] = 0.0;
1573: }
1574: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1575: nbcs = 0;
1576: if ((si+nx) == (M)) nbcs = ny;
1578: if (b) {
1579: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1580: VecAssemblyBegin(b);
1581: VecAssemblyEnd(b);
1582: }
1583: if (A) {
1584: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1585: }
1587: PetscFree(bc_vals);
1588: PetscFree(bc_global_ids);
1590: DMDAVecRestoreArray(cda,coords,&_coords);
1591: return(0);
1592: }
1594: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
1595: {
1596: DM cda;
1597: Vec coords;
1598: PetscInt si,sj,nx,ny,i,j;
1599: PetscInt M,N;
1600: DMDACoor2d **_coords;
1601: const PetscInt *g_idx;
1602: PetscInt *bc_global_ids;
1603: PetscScalar *bc_vals;
1604: PetscInt nbcs;
1605: PetscInt n_dofs;
1606: PetscErrorCode ierr;
1607: ISLocalToGlobalMapping ltogm;
1610: DMGetLocalToGlobalMapping(da,<ogm);
1611: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1613: DMGetCoordinateDM(da,&cda);
1614: DMGetCoordinatesLocal(da,&coords);
1615: DMDAVecGetArray(cda,coords,&_coords);
1616: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1617: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1619: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1620: PetscMalloc1(ny*n_dofs,&bc_vals);
1622: /* init the entries to -1 so VecSetValues will ignore them */
1623: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1625: i = 0;
1626: for (j = 0; j < ny; j++) {
1627: PetscInt local_id;
1629: local_id = i+j*nx;
1631: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1633: bc_vals[j] = 0.0;
1634: }
1635: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1636: nbcs = 0;
1637: if (si == 0) nbcs = ny;
1639: if (b) {
1640: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1641: VecAssemblyBegin(b);
1642: VecAssemblyEnd(b);
1643: }
1645: if (A) {
1646: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1647: }
1649: PetscFree(bc_vals);
1650: PetscFree(bc_global_ids);
1652: DMDAVecRestoreArray(cda,coords,&_coords);
1653: return(0);
1654: }
1656: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
1657: {
1658: DM cda;
1659: Vec coords;
1660: PetscInt si,sj,nx,ny,i,j;
1661: PetscInt M,N;
1662: DMDACoor2d **_coords;
1663: const PetscInt *g_idx;
1664: PetscInt *bc_global_ids;
1665: PetscScalar *bc_vals;
1666: PetscInt nbcs;
1667: PetscInt n_dofs;
1668: PetscErrorCode ierr;
1669: ISLocalToGlobalMapping ltogm;
1672: DMGetLocalToGlobalMapping(da,<ogm);
1673: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1675: DMGetCoordinateDM(da,&cda);
1676: DMGetCoordinatesLocal(da,&coords);
1677: DMDAVecGetArray(cda,coords,&_coords);
1678: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1679: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1681: PetscMalloc1(nx,&bc_global_ids);
1682: PetscMalloc1(nx,&bc_vals);
1684: /* init the entries to -1 so VecSetValues will ignore them */
1685: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1687: j = ny-1;
1688: for (i = 0; i < nx; i++) {
1689: PetscInt local_id;
1691: local_id = i+j*nx;
1693: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1695: bc_vals[i] = 0.0;
1696: }
1697: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1698: nbcs = 0;
1699: if ((sj+ny) == (N)) nbcs = nx;
1701: if (b) {
1702: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1703: VecAssemblyBegin(b);
1704: VecAssemblyEnd(b);
1705: }
1706: if (A) {
1707: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1708: }
1710: PetscFree(bc_vals);
1711: PetscFree(bc_global_ids);
1713: DMDAVecRestoreArray(cda,coords,&_coords);
1714: return(0);
1715: }
1717: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
1718: {
1719: DM cda;
1720: Vec coords;
1721: PetscInt si,sj,nx,ny,i,j;
1722: PetscInt M,N;
1723: DMDACoor2d **_coords;
1724: const PetscInt *g_idx;
1725: PetscInt *bc_global_ids;
1726: PetscScalar *bc_vals;
1727: PetscInt nbcs;
1728: PetscInt n_dofs;
1729: PetscErrorCode ierr;
1730: ISLocalToGlobalMapping ltogm;
1733: DMGetLocalToGlobalMapping(da,<ogm);
1734: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1736: DMGetCoordinateDM(da,&cda);
1737: DMGetCoordinatesLocal(da,&coords);
1738: DMDAVecGetArray(cda,coords,&_coords);
1739: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1740: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1742: PetscMalloc1(nx,&bc_global_ids);
1743: PetscMalloc1(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 = 0;
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] = 0.0;
1757: }
1758: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1759: nbcs = 0;
1760: if (sj == 0) nbcs = nx;
1762: if (b) {
1763: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1764: VecAssemblyBegin(b);
1765: VecAssemblyEnd(b);
1766: }
1767: if (A) {
1768: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1769: }
1771: PetscFree(bc_vals);
1772: PetscFree(bc_global_ids);
1774: DMDAVecRestoreArray(cda,coords,&_coords);
1775: return(0);
1776: }
1778: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
1779: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1780: {
1784: BCApplyZero_NORTH(da_Stokes,1,A,f);
1785: BCApplyZero_EAST(da_Stokes,0,A,f);
1786: BCApplyZero_SOUTH(da_Stokes,1,A,f);
1787: BCApplyZero_WEST(da_Stokes,0,A,f);
1788: return(0);
1789: }