Actual source code: ex43.c
1: static char help[] =
2: "Solves the incompressible, variable viscosity stokes equation in 2d on the unit domain \n\
3: using Q1Q1 elements, stabilized with Bochev's polynomial projection method. \n\
4: The models defined utilise free slip boundary conditions on all sides. \n\
5: Options: \n\
6: -mx : number elements in x-direciton \n\
7: -my : number elements in y-direciton \n\
8: -c_str : indicates the structure of the coefficients to use. \n\
9: -c_str 0 => Setup for an analytic solution with a vertical jump in viscosity. This problem is driven by the \n\
10: forcing function f = ( 0, sin(n_z pi y )cos( pi x ). \n\
11: Parameters: \n\
12: -solcx_eta0 : the viscosity to the left of the interface \n\
13: -solcx_eta1 : the viscosity to the right of the interface \n\
14: -solcx_xc : the location of the interface \n\
15: -solcx_nz : the wavenumber in the y direction \n\
16: -c_str 1 => Setup for a rectangular blob located in the origin of the domain. \n\
17: Parameters: \n\
18: -sinker_eta0 : the viscosity of the background fluid \n\
19: -sinker_eta1 : the viscosity of the blob \n\
20: -sinker_dx : the width of the blob \n\
21: -sinker_dy : the width of the blob \n\
22: -c_str 2 => Setup for a circular blob located in the origin of the domain. \n\
23: Parameters: \n\
24: -sinker_eta0 : the viscosity of the background fluid \n\
25: -sinker_eta1 : the viscosity of the blob \n\
26: -sinker_r : radius of the blob \n\
27: -use_gp_coords : evaluate the viscosity and the body force at the global coordinates of the quadrature points.\n\
28: By default, eta and the body force are evaulated at the element center and applied as a constant over the entire element.\n";
30: /* Contributed by Dave May */
32: #include <petscksp.h>
33: #include <petscdmda.h>
35: /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
36: #include "ex43-solCx.h"
38: static PetscErrorCode DMDABCApplyFreeSlip(DM,Mat,Vec);
41: #define NSD 2 /* number of spatial dimensions */
42: #define NODES_PER_EL 4 /* nodes per element */
43: #define U_DOFS 2 /* degrees of freedom per velocity node */
44: #define P_DOFS 1 /* degrees of freedom per pressure node */
45: #define GAUSS_POINTS 4
47: /* cell based evaluation */
48: typedef struct {
49: PetscScalar eta,fx,fy;
50: } Coefficients;
52: /* Gauss point based evaluation 8+4+4+4 = 20 */
53: typedef struct {
54: PetscScalar gp_coords[2*GAUSS_POINTS];
55: PetscScalar eta[GAUSS_POINTS];
56: PetscScalar fx[GAUSS_POINTS];
57: PetscScalar fy[GAUSS_POINTS];
58: } GaussPointCoefficients;
60: typedef struct {
61: PetscScalar u_dof;
62: PetscScalar v_dof;
63: PetscScalar p_dof;
64: } StokesDOF;
67: /*
69: D = [ 2.eta 0 0 ]
70: [ 0 2.eta 0 ]
71: [ 0 0 eta ]
73: B = [ d_dx 0 ]
74: [ 0 d_dy ]
75: [ d_dy d_dx ]
77: */
79: /* FEM routines */
80: /*
81: Element: Local basis function ordering
82: 1-----2
83: | |
84: | |
85: 0-----3
86: */
87: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
88: {
89: PetscScalar xi = _xi[0];
90: PetscScalar eta = _xi[1];
92: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
93: Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
94: Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
95: Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
96: }
98: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
99: {
100: PetscScalar xi = _xi[0];
101: PetscScalar eta = _xi[1];
103: GNi[0][0] = -0.25*(1.0-eta);
104: GNi[0][1] = -0.25*(1.0+eta);
105: GNi[0][2] = 0.25*(1.0+eta);
106: GNi[0][3] = 0.25*(1.0-eta);
108: GNi[1][0] = -0.25*(1.0-xi);
109: GNi[1][1] = 0.25*(1.0-xi);
110: GNi[1][2] = 0.25*(1.0+xi);
111: GNi[1][3] = -0.25*(1.0+xi);
112: }
114: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
115: {
116: PetscScalar J00,J01,J10,J11,J;
117: PetscScalar iJ00,iJ01,iJ10,iJ11;
118: PetscInt i;
120: J00 = J01 = J10 = J11 = 0.0;
121: for (i = 0; i < NODES_PER_EL; i++) {
122: PetscScalar cx = coords[ 2*i+0 ];
123: PetscScalar cy = coords[ 2*i+1 ];
125: J00 = J00+GNi[0][i]*cx; /* J_xx = dx/dxi */
126: J01 = J01+GNi[0][i]*cy; /* J_xy = dy/dxi */
127: J10 = J10+GNi[1][i]*cx; /* J_yx = dx/deta */
128: J11 = J11+GNi[1][i]*cy; /* J_yy = dy/deta */
129: }
130: J = (J00*J11)-(J01*J10);
132: iJ00 = J11/J;
133: iJ01 = -J01/J;
134: iJ10 = -J10/J;
135: iJ11 = J00/J;
138: for (i = 0; i < NODES_PER_EL; i++) {
139: GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
140: GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
141: }
143: if (det_J != NULL) {
144: *det_J = J;
145: }
146: }
148: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
149: {
150: *ngp = 4;
151: gp_xi[0][0] = -0.57735026919;gp_xi[0][1] = -0.57735026919;
152: gp_xi[1][0] = -0.57735026919;gp_xi[1][1] = 0.57735026919;
153: gp_xi[2][0] = 0.57735026919;gp_xi[2][1] = 0.57735026919;
154: gp_xi[3][0] = 0.57735026919;gp_xi[3][1] = -0.57735026919;
155: gp_weight[0] = 1.0;
156: gp_weight[1] = 1.0;
157: gp_weight[2] = 1.0;
158: gp_weight[3] = 1.0;
159: }
162: /* procs to the left claim the ghost node as their element */
165: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
166: {
167: PetscInt m,n,p,M,N,P;
168: PetscInt sx,sy,sz;
171: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
172: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
174: if (mxl != PETSC_NULL) {
175: *mxl = m;
176: if ((sx+m) == M) { /* last proc */
177: *mxl = m-1;
178: }
179: }
180: if (myl != PETSC_NULL) {
181: *myl = n;
182: if ((sy+n) == N) { /* last proc */
183: *myl = n-1;
184: }
185: }
186: if (mzl != PETSC_NULL) {
187: *mzl = p;
188: if ((sz+p) == P) { /* last proc */
189: *mzl = p-1;
190: }
191: }
192: return(0);
193: }
197: static PetscErrorCode DMDAGetElementCorners(DM da,
198: PetscInt *sx,PetscInt *sy,PetscInt *sz,
199: PetscInt *mx,PetscInt *my,PetscInt *mz)
200: {
201: PetscInt si,sj,sk;
204: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
206: *sx = si;
207: if (si != 0) {
208: *sx = si+1;
209: }
211: *sy = sj;
212: if (sj != 0) {
213: *sy = sj+1;
214: }
216: if (sk != PETSC_NULL) {
217: *sz = sk;
218: if (sk != 0) {
219: *sz = sk+1;
220: }
221: }
223: DMDAGetLocalElementSize(da,mx,my,mz);
225: return(0);
226: }
228: /*
229: i,j are the element indices
230: The unknown is a vector quantity.
231: The s[].c is used to indicate the degree of freedom.
232: */
235: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
236: {
238: /* velocity */
239: /* node 0 */
240: s_u[0].i = i;s_u[0].j = j;s_u[0].c = 0; /* Vx0 */
241: s_u[1].i = i;s_u[1].j = j;s_u[1].c = 1; /* Vy0 */
243: /* node 1 */
244: s_u[2].i = i;s_u[2].j = j+1;s_u[2].c = 0; /* Vx1 */
245: s_u[3].i = i;s_u[3].j = j+1;s_u[3].c = 1; /* Vy1 */
247: /* node 2 */
248: s_u[4].i = i+1;s_u[4].j = j+1;s_u[4].c = 0; /* Vx2 */
249: s_u[5].i = i+1;s_u[5].j = j+1;s_u[5].c = 1; /* Vy2 */
251: /* node 3 */
252: s_u[6].i = i+1;s_u[6].j = j;s_u[6].c = 0; /* Vx3 */
253: s_u[7].i = i+1;s_u[7].j = j;s_u[7].c = 1; /* Vy3 */
256: /* pressure */
257: s_p[0].i = i;s_p[0].j = j;s_p[0].c = 2; /* P0 */
258: s_p[1].i = i;s_p[1].j = j+1;s_p[1].c = 2; /* P0 */
259: s_p[2].i = i+1;s_p[2].j = j+1;s_p[2].c = 2; /* P1 */
260: s_p[3].i = i+1;s_p[3].j = j;s_p[3].c = 2; /* P1 */
261: return(0);
262: }
266: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
267: {
269: PetscMPIInt rank;
270: PetscInt proc_I,proc_J;
271: PetscInt cpu_x,cpu_y;
272: PetscInt local_mx,local_my;
273: Vec vlx,vly;
274: PetscInt *LX,*LY,i;
275: PetscScalar *_a;
276: Vec V_SEQ;
277: VecScatter ctx;
280: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
282: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
284: proc_J = rank/cpu_x;
285: proc_I = rank-cpu_x*proc_J;
287: PetscMalloc(sizeof(PetscInt)*cpu_x,&LX);
288: PetscMalloc(sizeof(PetscInt)*cpu_y,&LY);
290: DMDAGetLocalElementSize(da,&local_mx,&local_my,PETSC_NULL);
291: VecCreate(PETSC_COMM_WORLD,&vlx);
292: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
293: VecSetFromOptions(vlx);
295: VecCreate(PETSC_COMM_WORLD,&vly);
296: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
297: VecSetFromOptions(vly);
299: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
300: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
301: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
302: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
306: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
307: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
308: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
309: VecGetArray(V_SEQ,&_a);
310: for (i = 0; i < cpu_x; i++) {
311: LX[i] = (PetscInt)PetscRealPart(_a[i]);
312: }
313: VecRestoreArray(V_SEQ,&_a);
314: VecScatterDestroy(&ctx);
315: VecDestroy(&V_SEQ);
317: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
318: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
319: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
320: VecGetArray(V_SEQ,&_a);
321: for (i = 0; i < cpu_y; i++) {
322: LY[i] = (PetscInt)PetscRealPart(_a[i]);
323: }
324: VecRestoreArray(V_SEQ,&_a);
325: VecScatterDestroy(&ctx);
326: VecDestroy(&V_SEQ);
330: *_lx = LX;
331: *_ly = LY;
333: VecDestroy(&vlx);
334: VecDestroy(&vly);
336: return(0);
337: }
341: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
342: {
343: DM cda;
344: Vec coords;
345: DMDACoor2d **_coords;
346: PetscInt si,sj,nx,ny,i,j;
347: FILE *fp;
348: char fname[PETSC_MAX_PATH_LEN];
349: PetscMPIInt rank;
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,"### Element geometry for processor %1.4d ### \n",rank);
360: DMDAGetCoordinateDA(da,&cda);
361: DMDAGetGhostedCoordinates(da,&coords);
362: DMDAVecGetArray(cda,coords,&_coords);
363: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
364: for (j = sj; j < sj+ny-1; j++) {
365: for (i = si; i < si+nx-1; i++) {
366: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
367: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i].x),PetscRealPart(_coords[j+1][i].y));
368: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j+1][i+1].x),PetscRealPart(_coords[j+1][i+1].y));
369: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",PetscRealPart(_coords[j][i+1].x),PetscRealPart(_coords[j][i+1].y));
370: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",PetscRealPart(_coords[j][i].x),PetscRealPart(_coords[j][i].y));
371: }
372: }
373: DMDAVecRestoreArray(cda,coords,&_coords);
375: PetscFClose(PETSC_COMM_SELF,fp);
376: return(0);
377: }
381: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
382: {
383: DM cda;
384: Vec coords,local_fields;
385: DMDACoor2d **_coords;
386: FILE *fp;
387: char fname[PETSC_MAX_PATH_LEN];
388: PetscMPIInt rank;
389: PetscInt si,sj,nx,ny,i,j;
390: PetscInt n_dofs,d;
391: PetscScalar *_fields;
395: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
396: PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
397: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
398: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
400: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
401: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
402: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
403: for (d = 0; d < n_dofs; d++) {
404: const char *field_name;
405: DMDAGetFieldName(da,d,&field_name);
406: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
407: }
408: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
411: DMDAGetCoordinateDA(da,&cda);
412: DMDAGetGhostedCoordinates(da,&coords);
413: DMDAVecGetArray(cda,coords,&_coords);
414: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
416: DMCreateLocalVector(da,&local_fields);
417: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
418: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
419: VecGetArray(local_fields,&_fields);
422: for (j = sj; j < sj+ny; j++) {
423: for (i = si; i < si+nx; i++) {
424: PetscScalar coord_x,coord_y;
425: PetscScalar field_d;
427: coord_x = _coords[j][i].x;
428: coord_y = _coords[j][i].y;
430: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
431: for (d = 0; d < n_dofs; d++) {
432: field_d = _fields[ n_dofs*((i-si)+(j-sj)*(nx))+d ];
433: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",PetscRealPart(field_d));
434: }
435: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
436: }
437: }
438: VecRestoreArray(local_fields,&_fields);
439: VecDestroy(&local_fields);
441: DMDAVecRestoreArray(cda,coords,&_coords);
443: PetscFClose(PETSC_COMM_SELF,fp);
444: return(0);
445: }
449: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
450: {
451: DM cda;
452: Vec local_fields;
453: FILE *fp;
454: char fname[PETSC_MAX_PATH_LEN];
455: PetscMPIInt rank;
456: PetscInt si,sj,nx,ny,i,j,p;
457: PetscInt n_dofs,d;
458: GaussPointCoefficients **_coefficients;
459: PetscErrorCode ierr;
462: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
463: PetscSNPrintf(fname,sizeof fname,"%s-p%1.4d.dat",prefix,rank);
464: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
465: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
467: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
468: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
469: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
470: for (d = 0; d < n_dofs; d++) {
471: const char *field_name;
472: DMDAGetFieldName(da,d,&field_name);
473: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
474: }
475: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
478: DMDAGetCoordinateDA(da,&cda);
479: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
481: DMCreateLocalVector(da,&local_fields);
482: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
483: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
484: DMDAVecGetArray(da,local_fields,&_coefficients);
487: for (j = sj; j < sj+ny; j++) {
488: for (i = si; i < si+nx; i++) {
489: PetscScalar coord_x,coord_y;
491: for (p = 0; p < GAUSS_POINTS; p++) {
492: coord_x = _coefficients[j][i].gp_coords[2*p];
493: coord_y = _coefficients[j][i].gp_coords[2*p+1];
495: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",PetscRealPart(coord_x),PetscRealPart(coord_y));
497: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e %1.6e",PetscRealPart(_coefficients[j][i].eta[p]),PetscRealPart(_coefficients[j][i].fx[p]),PetscRealPart(_coefficients[j][i].fy[p]));
498: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
499: }
500: }
501: }
502: DMDAVecRestoreArray(da,local_fields,&_coefficients);
503: VecDestroy(&local_fields);
505: PetscFClose(PETSC_COMM_SELF,fp);
506: return(0);
507: }
510: static PetscInt ASS_MAP_wIwDI_uJuDJ(
511: PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,
512: PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
513: {
514: PetscInt ij;
515: PetscInt r,c,nc;
517: nc = u_NPE*u_dof;
519: r = w_dof*wi+wd;
520: c = u_dof*ui+ud;
522: ij = r*nc+c;
524: return ij;
525: }
527: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
528: {
529: PetscInt ngp;
530: PetscScalar gp_xi[GAUSS_POINTS][2];
531: PetscScalar gp_weight[GAUSS_POINTS];
532: PetscInt p,i,j,k;
533: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
534: PetscScalar J_p,tildeD[3];
535: PetscScalar B[3][U_DOFS*NODES_PER_EL];
538: /* define quadrature rule */
539: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
541: /* evaluate integral */
542: for (p = 0; p < ngp; p++) {
543: ConstructQ12D_GNi(gp_xi[p],GNi_p);
544: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
546: for (i = 0; i < NODES_PER_EL; i++) {
547: PetscScalar d_dx_i = GNx_p[0][i];
548: PetscScalar d_dy_i = GNx_p[1][i];
550: B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
551: B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
552: B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
553: }
556: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
557: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
558: tildeD[2] = gp_weight[p]*J_p*eta[p];
560: /* form Bt tildeD B */
561: /*
562: Ke_ij = Bt_ik . D_kl . B_lj
563: = B_ki . D_kl . B_lj
564: = B_ki . D_kk . B_kj
565: */
566: for (i = 0; i < 8; i++) {
567: for (j = 0; j < 8; j++) {
568: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
569: Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
570: }
571: }
572: }
573: }
574: }
576: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
577: {
578: PetscInt ngp;
579: PetscScalar gp_xi[GAUSS_POINTS][2];
580: PetscScalar gp_weight[GAUSS_POINTS];
581: PetscInt p,i,j,di;
582: PetscScalar Ni_p[NODES_PER_EL];
583: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
584: PetscScalar J_p,fac;
587: /* define quadrature rule */
588: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
590: /* evaluate integral */
591: for (p = 0; p < ngp; p++) {
592: ConstructQ12D_Ni(gp_xi[p],Ni_p);
593: ConstructQ12D_GNi(gp_xi[p],GNi_p);
594: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
595: fac = gp_weight[p]*J_p;
597: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
598: for (di = 0; di < NSD; di++) { /* u dofs */
599: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
600: PetscInt IJ;
601: /* Ke[4*u_idx+j] = Ke[4*u_idx+j] - GNx_p[di][i] * Ni_p[j] * fac; */
602: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);
604: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
605: }
606: }
607: }
608: }
609: }
611: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
612: {
613: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
614: PetscInt i,j;
615: PetscInt nr_g,nc_g;
617: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
618: FormGradientOperatorQ1(Ge,coords);
620: nr_g = U_DOFS*NODES_PER_EL;
621: nc_g = P_DOFS*NODES_PER_EL;
623: for (i = 0; i < nr_g; i++) {
624: for (j = 0; j < nc_g; j++) {
625: De[nr_g*j+i] = Ge[nc_g*i+j];
626: }
627: }
628: }
630: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
631: {
632: PetscInt ngp;
633: PetscScalar gp_xi[GAUSS_POINTS][2];
634: PetscScalar gp_weight[GAUSS_POINTS];
635: PetscInt p,i,j;
636: PetscScalar Ni_p[NODES_PER_EL];
637: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
638: PetscScalar J_p,fac,eta_avg;
641: /* define quadrature rule */
642: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
644: /* evaluate integral */
645: for (p = 0; p < ngp; p++) {
646: ConstructQ12D_Ni(gp_xi[p],Ni_p);
647: ConstructQ12D_GNi(gp_xi[p],GNi_p);
648: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
649: fac = gp_weight[p]*J_p;
651: for (i = 0; i < NODES_PER_EL; i++) {
652: for (j = 0; j < NODES_PER_EL; j++) {
653: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
654: }
655: }
656: }
658: /* scale */
659: eta_avg = 0.0;
660: for (p = 0; p < ngp; p++) {
661: eta_avg += eta[p];
662: }
663: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
664: fac = 1.0/eta_avg;
665: for (i = 0; i < NODES_PER_EL; i++) {
666: for (j = 0; j < NODES_PER_EL; j++) {
667: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
668: }
669: }
670: }
672: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
673: {
674: PetscInt ngp;
675: PetscScalar gp_xi[GAUSS_POINTS][2];
676: PetscScalar gp_weight[GAUSS_POINTS];
677: PetscInt p,i,j;
678: PetscScalar Ni_p[NODES_PER_EL];
679: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
680: PetscScalar J_p,fac,eta_avg;
683: /* define quadrature rule */
684: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
686: /* evaluate integral */
687: for (p = 0; p < ngp; p++) {
688: ConstructQ12D_Ni(gp_xi[p],Ni_p);
689: ConstructQ12D_GNi(gp_xi[p],GNi_p);
690: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
691: fac = gp_weight[p]*J_p;
693: for (i = 0; i < NODES_PER_EL; i++) {
694: for (j = 0; j < NODES_PER_EL; j++) {
695: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
696: }
697: }
698: }
700: /* scale */
701: eta_avg = 0.0;
702: for (p = 0; p < ngp; p++) {
703: eta_avg += eta[p];
704: }
705: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
706: fac = 1.0/eta_avg;
707: for (i = 0; i < NODES_PER_EL; i++) {
708: for (j = 0; j < NODES_PER_EL; j++) {
709: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
710: }
711: }
712: }
714: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
715: {
716: PetscInt ngp;
717: PetscScalar gp_xi[GAUSS_POINTS][2];
718: PetscScalar gp_weight[GAUSS_POINTS];
719: PetscInt p,i;
720: PetscScalar Ni_p[NODES_PER_EL];
721: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
722: PetscScalar J_p,fac;
725: /* define quadrature rule */
726: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
728: /* evaluate integral */
729: for (p = 0; p < ngp; p++) {
730: ConstructQ12D_Ni(gp_xi[p],Ni_p);
731: ConstructQ12D_GNi(gp_xi[p],GNi_p);
732: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
733: fac = gp_weight[p]*J_p;
735: for (i = 0; i < NODES_PER_EL; i++) {
736: Fe[NSD*i ] += fac*Ni_p[i]*fx[p];
737: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
738: }
739: }
740: }
744: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
745: {
747: /* get coords for the element */
748: el_coords[NSD*0+0] = _coords[ej][ei].x;el_coords[NSD*0+1] = _coords[ej][ei].y;
749: el_coords[NSD*1+0] = _coords[ej+1][ei].x;el_coords[NSD*1+1] = _coords[ej+1][ei].y;
750: el_coords[NSD*2+0] = _coords[ej+1][ei+1].x;el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
751: el_coords[NSD*3+0] = _coords[ej][ei+1].x;el_coords[NSD*3+1] = _coords[ej][ei+1].y;
752: return(0);
753: }
757: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
758: {
759: DM cda;
760: Vec coords;
761: DMDACoor2d **_coords;
762: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
763: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
764: PetscInt sex,sey,mx,my;
765: PetscInt ei,ej;
766: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
767: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
768: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
769: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
770: PetscScalar el_coords[NODES_PER_EL*NSD];
771: Vec local_properties;
772: GaussPointCoefficients **props;
773: PetscScalar *prop_eta;
774: PetscErrorCode ierr;
778: /* setup for coords */
779: DMDAGetCoordinateDA(stokes_da,&cda);
780: DMDAGetGhostedCoordinates(stokes_da,&coords);
781: DMDAVecGetArray(cda,coords,&_coords);
783: /* setup for coefficients */
784: DMCreateLocalVector(properties_da,&local_properties);
785: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
786: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
787: DMDAVecGetArray(properties_da,local_properties,&props);
789: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
790: for (ej = sey; ej < sey+my; ej++) {
791: for (ei = sex; ei < sex+mx; ei++) {
792: /* get coords for the element */
793: GetElementCoords(_coords,ei,ej,el_coords);
795: /* get coefficients for the element */
796: prop_eta = props[ej][ei].eta;
798: /* initialise element stiffness matrix */
799: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
800: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
801: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
802: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
804: /* form element stiffness matrix */
805: FormStressOperatorQ1(Ae,el_coords,prop_eta);
806: FormGradientOperatorQ1(Ge,el_coords);
807: FormDivergenceOperatorQ1(De,el_coords);
808: FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);
810: /* insert element matrix into global matrix */
811: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
812: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
813: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
814: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
815: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
816: }
817: }
818: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
819: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
821: DMDAVecRestoreArray(cda,coords,&_coords);
823: DMDAVecRestoreArray(properties_da,local_properties,&props);
824: VecDestroy(&local_properties);
825: return(0);
826: }
830: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
831: {
832: DM cda;
833: Vec coords;
834: DMDACoor2d **_coords;
835: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
836: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
837: PetscInt sex,sey,mx,my;
838: PetscInt ei,ej;
839: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
840: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
841: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
842: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
843: PetscScalar el_coords[NODES_PER_EL*NSD];
844: Vec local_properties;
845: GaussPointCoefficients **props;
846: PetscScalar *prop_eta;
847: PetscErrorCode ierr;
850: /* setup for coords */
851: DMDAGetCoordinateDA(stokes_da,&cda);
852: DMDAGetGhostedCoordinates(stokes_da,&coords);
853: DMDAVecGetArray(cda,coords,&_coords);
855: /* setup for coefficients */
856: DMCreateLocalVector(properties_da,&local_properties);
857: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
858: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
859: DMDAVecGetArray(properties_da,local_properties,&props);
861: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
862: for (ej = sey; ej < sey+my; ej++) {
863: for (ei = sex; ei < sex+mx; ei++) {
864: /* get coords for the element */
865: GetElementCoords(_coords,ei,ej,el_coords);
867: /* get coefficients for the element */
868: prop_eta = props[ej][ei].eta;
870: /* initialise element stiffness matrix */
871: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
872: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
873: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
874: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
877: /* form element stiffness matrix */
878: FormStressOperatorQ1(Ae,el_coords,prop_eta);
879: FormGradientOperatorQ1(Ge,el_coords);
880: /* FormDivergenceOperatorQ1( De, el_coords ); */
881: FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);
883: /* insert element matrix into global matrix */
884: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
885: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
886: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
887: /* MatSetValuesStencil( A, NODES_PER_EL*P_DOFS,p_eqn, NODES_PER_EL*U_DOFS,u_eqn, De, ADD_VALUES ); */
888: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
889: }
890: }
891: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
892: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
894: DMDAVecRestoreArray(cda,coords,&_coords);
896: DMDAVecRestoreArray(properties_da,local_properties,&props);
897: VecDestroy(&local_properties);
898: return(0);
899: }
903: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
904: {
905: PetscInt n;
908: for (n = 0; n < 4; n++) {
909: 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 ];
910: 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];
911: 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 ];
912: }
913: return(0);
914: }
918: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
919: {
920: DM cda;
921: Vec coords;
922: DMDACoor2d **_coords;
923: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
924: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
925: PetscInt sex,sey,mx,my;
926: PetscInt ei,ej;
927: PetscScalar Fe[NODES_PER_EL*U_DOFS];
928: PetscScalar He[NODES_PER_EL*P_DOFS];
929: PetscScalar el_coords[NODES_PER_EL*NSD];
930: Vec local_properties;
931: GaussPointCoefficients **props;
932: PetscScalar *prop_fx,*prop_fy;
933: Vec local_F;
934: StokesDOF **ff;
935: PetscErrorCode ierr;
938: /* setup for coords */
939: DMDAGetCoordinateDA(stokes_da,&cda);
940: DMDAGetGhostedCoordinates(stokes_da,&coords);
941: DMDAVecGetArray(cda,coords,&_coords);
943: /* setup for coefficients */
944: DMGetLocalVector(properties_da,&local_properties);
945: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
946: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
947: DMDAVecGetArray(properties_da,local_properties,&props);
949: /* get acces to the vector */
950: DMGetLocalVector(stokes_da,&local_F);
951: VecZeroEntries(local_F);
952: DMDAVecGetArray(stokes_da,local_F,&ff);
955: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
956: for (ej = sey; ej < sey+my; ej++) {
957: for (ei = sex; ei < sex+mx; ei++) {
958: /* get coords for the element */
959: GetElementCoords(_coords,ei,ej,el_coords);
961: /* get coefficients for the element */
962: prop_fx = props[ej][ei].fx;
963: prop_fy = props[ej][ei].fy;
965: /* initialise element stiffness matrix */
966: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
967: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
970: /* form element stiffness matrix */
971: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
973: /* insert element matrix into global matrix */
974: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
976: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
977: }
978: }
980: DMDAVecRestoreArray(stokes_da,local_F,&ff);
981: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
982: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
983: DMRestoreLocalVector(stokes_da,&local_F);
986: DMDAVecRestoreArray(cda,coords,&_coords);
988: DMDAVecRestoreArray(properties_da,local_properties,&props);
989: DMRestoreLocalVector(properties_da,&local_properties);
990: return(0);
991: }
995: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,
996: PetscInt mx,PetscInt my,
997: DM *_da,Vec *_X)
998: {
999: DM da,cda;
1000: Vec X,local_X;
1001: StokesDOF **_stokes;
1002: Vec coords;
1003: DMDACoor2d **_coords;
1004: PetscInt si,sj,ei,ej,i,j;
1008: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1009: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,3,1,PETSC_NULL,PETSC_NULL,&da);
1010: DMDASetFieldName(da,0,"anlytic_Vx");
1011: DMDASetFieldName(da,1,"anlytic_Vy");
1012: DMDASetFieldName(da,2,"analytic_P");
1015: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,PETSC_NULL,PETSC_NULL);
1018: DMDAGetGhostedCoordinates(da,&coords);
1019: DMDAGetCoordinateDA(da,&cda);
1020: DMDAVecGetArray(cda,coords,&_coords);
1022: DMCreateGlobalVector(da,&X);
1023: DMCreateLocalVector(da,&local_X);
1024: DMDAVecGetArray(da,local_X,&_stokes);
1026: DMDAGetGhostCorners(da,&si,&sj,0,&ei,&ej,0);
1027: for (j = sj; j < sj+ej; j++) {
1028: for (i = si; i < si+ei; i++) {
1029: double pos[2],pressure,vel[2],total_stress[3],strain_rate[3];
1031: pos[0] = PetscRealPart(_coords[j][i].x);
1032: pos[1] = PetscRealPart(_coords[j][i].y);
1034: evaluate_solCx(pos,eta0,eta1,xc,nz,vel,&pressure,total_stress,strain_rate);
1036: _stokes[j][i].u_dof = vel[0];
1037: _stokes[j][i].v_dof = vel[1];
1038: _stokes[j][i].p_dof = pressure;
1039: }
1040: }
1041: DMDAVecRestoreArray(da,local_X,&_stokes);
1042: DMDAVecRestoreArray(cda,coords,&_coords);
1044: DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1045: DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);
1047: VecDestroy(&local_X);
1049: *_da = da;
1050: *_X = X;
1052: return(0);
1053: }
1057: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
1058: {
1060: /* get the nodal fields */
1061: 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;
1062: 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;
1063: 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;
1064: 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;
1065: return(0);
1066: }
1070: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
1071: {
1072: DM cda;
1073: Vec coords,X_analytic_local,X_local;
1074: DMDACoor2d **_coords;
1075: PetscInt sex,sey,mx,my;
1076: PetscInt ei,ej;
1077: PetscScalar el_coords[NODES_PER_EL*NSD];
1078: StokesDOF **stokes_analytic,**stokes;
1079: StokesDOF stokes_analytic_e[4],stokes_e[4];
1081: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1082: PetscScalar Ni_p[NODES_PER_EL];
1083: PetscInt ngp;
1084: PetscScalar gp_xi[GAUSS_POINTS][2];
1085: PetscScalar gp_weight[GAUSS_POINTS];
1086: PetscInt p,i;
1087: PetscScalar J_p,fac;
1088: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1089: PetscInt M;
1090: PetscReal xymin[2],xymax[2];
1094: /* define quadrature rule */
1095: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1097: /* setup for coords */
1098: DMDAGetCoordinateDA(stokes_da,&cda);
1099: DMDAGetGhostedCoordinates(stokes_da,&coords);
1100: DMDAVecGetArray(cda,coords,&_coords);
1102: /* setup for analytic */
1103: DMCreateLocalVector(stokes_da,&X_analytic_local);
1104: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1105: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1106: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1108: /* setup for solution */
1109: DMCreateLocalVector(stokes_da,&X_local);
1110: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1111: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1112: DMDAVecGetArray(stokes_da,X_local,&stokes);
1114: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1115: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1117: h = (xymax[0]-xymin[0])/((double)M);
1119: tp_L2 = tu_L2 = tu_H1 = 0.0;
1121: DMDAGetElementCorners(stokes_da,&sex,&sey,0,&mx,&my,0);
1122: for (ej = sey; ej < sey+my; ej++) {
1123: for (ei = sex; ei < sex+mx; ei++) {
1124: /* get coords for the element */
1125: GetElementCoords(_coords,ei,ej,el_coords);
1126: StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1127: StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);
1129: /* evaluate integral */
1130: p_e_L2 = 0.0;
1131: u_e_L2 = 0.0;
1132: u_e_H1 = 0.0;
1133: for (p = 0; p < ngp; p++) {
1134: ConstructQ12D_Ni(gp_xi[p],Ni_p);
1135: ConstructQ12D_GNi(gp_xi[p],GNi_p);
1136: ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1137: fac = gp_weight[p]*J_p;
1139: for (i = 0; i < NODES_PER_EL; i++) {
1140: PetscScalar u_error,v_error;
1142: 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);
1144: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1145: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1146: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);
1148: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1149: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error /* du/dy */
1150: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1151: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error); /* dv/dy */
1152: }
1153: }
1155: tp_L2 += p_e_L2;
1156: tu_L2 += u_e_L2;
1157: tu_H1 += u_e_H1;
1158: }
1159: }
1160: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1161: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1162: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1163: p_L2 = PetscSqrtScalar(p_L2);
1164: u_L2 = PetscSqrtScalar(u_L2);
1165: u_H1 = PetscSqrtScalar(u_H1);
1167: PetscPrintf(PETSC_COMM_WORLD,"%1.4e %1.4e %1.4e %1.4e \n",PetscRealPart(h),PetscRealPart(p_L2),PetscRealPart(u_L2),PetscRealPart(u_H1));
1170: DMDAVecRestoreArray(cda,coords,&_coords);
1172: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1173: VecDestroy(&X_analytic_local);
1174: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1175: VecDestroy(&X_local);
1176: return(0);
1177: }
1181: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1182: {
1183: DM da_Stokes,da_prop;
1184: PetscInt u_dof,p_dof,dof,stencil_width;
1185: Mat A,B;
1186: PetscInt mxl,myl;
1187: DM prop_cda,vel_cda;
1188: Vec prop_coords,vel_coords;
1189: PetscInt si,sj,nx,ny,i,j,p;
1190: Vec f,X;
1191: PetscInt prop_dof,prop_stencil_width;
1192: Vec properties,l_properties;
1193: PetscReal dx,dy;
1194: PetscInt M,N;
1195: DMDACoor2d **_prop_coords,**_vel_coords;
1196: GaussPointCoefficients **element_props;
1197: PetscInt its;
1198: KSP ksp_S;
1199: PetscInt coefficient_structure = 0;
1200: PetscInt cpu_x,cpu_y,*lx = PETSC_NULL,*ly = PETSC_NULL;
1201: PetscBool use_gp_coords = PETSC_FALSE;
1202: PetscErrorCode ierr;
1205: /* Generate the da for velocity and pressure */
1206: /*
1207: We use Q1 elements for the temperature.
1208: FEM has a 9-point stencil (BOX) or connectivity pattern
1209: Num nodes in each direction is mx+1, my+1
1210: */
1211: u_dof = U_DOFS; /* Vx, Vy - velocities */
1212: p_dof = P_DOFS; /* p - pressure */
1213: dof = u_dof+p_dof;
1214: stencil_width = 1;
1215: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1216: mx+1,my+1,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,PETSC_NULL,PETSC_NULL,&da_Stokes);
1217: DMDASetFieldName(da_Stokes,0,"Vx");
1218: DMDASetFieldName(da_Stokes,1,"Vy");
1219: DMDASetFieldName(da_Stokes,2,"P");
1221: /* unit box [0,1] x [0,1] */
1222: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,PETSC_NULL,PETSC_NULL);
1225: /* Generate element properties, we will assume all material properties are constant over the element */
1226: /* local number of elements */
1227: DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,PETSC_NULL);
1229: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! // */
1230: DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1231: DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);
1233: prop_dof = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1234: prop_stencil_width = 0;
1235: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1236: mx,my,cpu_x,cpu_y,prop_dof,prop_stencil_width,lx,ly,&da_prop);
1237: PetscFree(lx);
1238: PetscFree(ly);
1240: /* define centroid positions */
1241: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1242: dx = 1.0/((PetscReal)(M));
1243: dy = 1.0/((PetscReal)(N));
1245: DMDASetUniformCoordinates(da_prop,0.0+0.5*dx,1.0-0.5*dx,0.0+0.5*dy,1.0-0.5*dy,PETSC_NULL,PETSC_NULL);
1247: /* define coefficients */
1248: PetscOptionsGetInt(PETSC_NULL,"-c_str",&coefficient_structure,PETSC_NULL);
1249: /* PetscPrintf( PETSC_COMM_WORLD, "Using coeficient structure %D \n", coefficient_structure ); */
1251: DMCreateGlobalVector(da_prop,&properties);
1252: DMCreateLocalVector(da_prop,&l_properties);
1253: DMDAVecGetArray(da_prop,l_properties,&element_props);
1255: DMDAGetCoordinateDA(da_prop,&prop_cda);
1256: DMDAGetGhostedCoordinates(da_prop,&prop_coords);
1257: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
1259: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
1261: DMDAGetCoordinateDA(da_Stokes,&vel_cda);
1262: DMDAGetGhostedCoordinates(da_Stokes,&vel_coords);
1263: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1266: /* interpolate the coordinates */
1267: for (j = sj; j < sj+ny; j++) {
1268: for (i = si; i < si+nx; i++) {
1269: PetscInt ngp;
1270: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1271: PetscScalar el_coords[8];
1273: GetElementCoords(_vel_coords,i,j,el_coords);
1274: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1276: for (p = 0; p < GAUSS_POINTS; p++) {
1277: PetscScalar gp_x,gp_y;
1278: PetscInt n;
1279: PetscScalar xi_p[2],Ni_p[4];
1281: xi_p[0] = gp_xi[p][0];
1282: xi_p[1] = gp_xi[p][1];
1283: ConstructQ12D_Ni(xi_p,Ni_p);
1285: gp_x = 0.0;
1286: gp_y = 0.0;
1287: for (n = 0; n < NODES_PER_EL; n++) {
1288: gp_x = gp_x+Ni_p[n]*el_coords[2*n ];
1289: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1290: }
1291: element_props[j][i].gp_coords[2*p ] = gp_x;
1292: element_props[j][i].gp_coords[2*p+1] = gp_y;
1293: }
1294: }
1295: }
1297: /* define the coefficients */
1298: PetscOptionsGetBool(PETSC_NULL,"-use_gp_coords",&use_gp_coords,0);
1300: for (j = sj; j < sj+ny; j++) {
1301: for (i = si; i < si+nx; i++) {
1302: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1303: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1304: PetscReal coord_x,coord_y;
1306: if (coefficient_structure == 0) {
1307: PetscReal opts_eta0,opts_eta1,opts_xc;
1308: PetscInt opts_nz;
1310: opts_eta0 = 1.0;
1311: opts_eta1 = 1.0;
1312: opts_xc = 0.5;
1313: opts_nz = 1;
1315: PetscOptionsGetReal(PETSC_NULL,"-solcx_eta0",&opts_eta0,0);
1316: PetscOptionsGetReal(PETSC_NULL,"-solcx_eta1",&opts_eta1,0);
1317: PetscOptionsGetReal(PETSC_NULL,"-solcx_xc",&opts_xc,0);
1318: PetscOptionsGetInt(PETSC_NULL,"-solcx_nz",&opts_nz,0);
1320: for (p = 0; p < GAUSS_POINTS; p++) {
1321: coord_x = centroid_x;
1322: coord_y = centroid_y;
1323: if (use_gp_coords == PETSC_TRUE) {
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: }
1329: element_props[j][i].eta[p] = opts_eta0;
1330: if (coord_x > opts_xc) {
1331: element_props[j][i].eta[p] = opts_eta1;
1332: }
1334: element_props[j][i].fx[p] = 0.0;
1335: element_props[j][i].fy[p] = sin((double)opts_nz*PETSC_PI*coord_y)*cos(1.0*PETSC_PI*coord_x);
1336: }
1337: } else if (coefficient_structure == 1) { /* square sinker */
1338: PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;
1340: opts_eta0 = 1.0;
1341: opts_eta1 = 1.0;
1342: opts_dx = 0.50;
1343: opts_dy = 0.50;
1345: PetscOptionsGetReal(PETSC_NULL,"-sinker_eta0",&opts_eta0,0);
1346: PetscOptionsGetReal(PETSC_NULL,"-sinker_eta1",&opts_eta1,0);
1347: PetscOptionsGetReal(PETSC_NULL,"-sinker_dx",&opts_dx,0);
1348: PetscOptionsGetReal(PETSC_NULL,"-sinker_dy",&opts_dy,0);
1351: for (p = 0; p < GAUSS_POINTS; p++) {
1352: coord_x = centroid_x;
1353: coord_y = centroid_y;
1354: if (use_gp_coords == PETSC_TRUE) {
1355: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1356: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1357: }
1359: element_props[j][i].eta[p] = opts_eta0;
1360: element_props[j][i].fx[p] = 0.0;
1361: element_props[j][i].fy[p] = 0.0;
1363: if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1364: if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1365: element_props[j][i].eta[p] = opts_eta1;
1366: element_props[j][i].fx[p] = 0.0;
1367: element_props[j][i].fy[p] = -1.0;
1368: }
1369: }
1370: }
1371: } else if (coefficient_structure == 2) { /* circular sinker */
1372: PetscReal opts_eta0,opts_eta1,opts_r,radius2;
1374: opts_eta0 = 1.0;
1375: opts_eta1 = 1.0;
1376: opts_r = 0.25;
1378: PetscOptionsGetReal(PETSC_NULL,"-sinker_eta0",&opts_eta0,0);
1379: PetscOptionsGetReal(PETSC_NULL,"-sinker_eta1",&opts_eta1,0);
1380: PetscOptionsGetReal(PETSC_NULL,"-sinker_r",&opts_r,0);
1382: for (p = 0; p < GAUSS_POINTS; p++) {
1383: coord_x = centroid_x;
1384: coord_y = centroid_y;
1385: if (use_gp_coords == PETSC_TRUE) {
1386: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1387: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1388: }
1390: element_props[j][i].eta[p] = opts_eta0;
1391: element_props[j][i].fx[p] = 0.0;
1392: element_props[j][i].fy[p] = 0.0;
1394: radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1395: if (radius2 < opts_r*opts_r) {
1396: element_props[j][i].eta[p] = opts_eta1;
1397: element_props[j][i].fx[p] = 0.0;
1398: element_props[j][i].fy[p] = -1.0;
1399: }
1400: }
1401: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1402: }
1403: }
1404: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
1406: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1408: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1409: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1410: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
1413: DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1414: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coeffcients for Stokes eqn.","properties");
1417: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1418: DMGetMatrix(da_Stokes,MATAIJ,&A);
1419: DMGetMatrix(da_Stokes,MATAIJ,&B);
1420: MatGetVecs(A,&f,&X);
1422: /* assemble A11 */
1423: MatZeroEntries(A);
1424: MatZeroEntries(B);
1425: VecZeroEntries(f);
1427: AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1428: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1429: /* build force vector */
1430: AssembleF_Stokes(f,da_Stokes,da_prop,properties);
1432: DMDABCApplyFreeSlip(da_Stokes,A,f);
1433: DMDABCApplyFreeSlip(da_Stokes,B,PETSC_NULL);
1435: /* SOLVE */
1436: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1437: KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
1438: KSPSetOperators(ksp_S,A,B,SAME_NONZERO_PATTERN);
1439: KSPSetFromOptions(ksp_S);
1440: {
1441: PC pc;
1442: const PetscInt ufields[] = {0,1},pfields[1] = {2};
1443: KSPGetPC(ksp_S,&pc);
1444: PCFieldSplitSetFields(pc,"u",2,ufields);
1445: PCFieldSplitSetFields(pc,"p",1,pfields);
1446: }
1448: KSPSolve(ksp_S,f,X);
1449: DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
1451: KSPGetIterationNumber(ksp_S,&its);
1453: /*
1454: {
1455: Vec x;
1456: PetscInt L;
1457: VecDuplicate(f,&x);
1458: MatMult(A,X, x );
1459: VecAXPY( x, -1.0, f );
1460: VecNorm( x, NORM_2, &nrm );
1461: PetscPrintf( PETSC_COMM_WORLD, "Its. %1.4d, norm |AX-f| = %1.5e \n", its, nrm );
1462: VecDestroy(&x);
1464: VecNorm( X, NORM_2, &nrm );
1465: VecGetSize( X, &L );
1466: PetscPrintf( PETSC_COMM_WORLD, " norm |X|/sqrt{N} = %1.5e \n", nrm/sqrt( (PetscScalar)L ) );
1467: }
1468: */
1470: if (coefficient_structure == 0) {
1471: PetscReal opts_eta0,opts_eta1,opts_xc;
1472: PetscInt opts_nz,N;
1473: DM da_Stokes_analytic;
1474: Vec X_analytic;
1475: PetscReal nrm1[3],nrm2[3],nrmI[3];
1477: opts_eta0 = 1.0;
1478: opts_eta1 = 1.0;
1479: opts_xc = 0.5;
1480: opts_nz = 1;
1482: PetscOptionsGetReal(PETSC_NULL,"-solcx_eta0",&opts_eta0,0);
1483: PetscOptionsGetReal(PETSC_NULL,"-solcx_eta1",&opts_eta1,0);
1484: PetscOptionsGetReal(PETSC_NULL,"-solcx_xc",&opts_xc,0);
1485: PetscOptionsGetInt(PETSC_NULL,"-solcx_nz",&opts_nz,0);
1488: DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1489: DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
1491: DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);
1494: VecAXPY(X_analytic,-1.0,X);
1495: VecGetSize(X_analytic,&N);
1496: N = N/3;
1498: VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1499: VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1500: VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);
1502: VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1503: VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1504: VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);
1506: VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1507: VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1508: VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);
1510: /*
1511: PetscPrintf( PETSC_COMM_WORLD, "%1.4e %1.4e %1.4e %1.4e %1.4e %1.4e %1.4e %1.4e %1.4e %1.4e\n",
1512: 1.0/mx,
1513: nrm1[0]/(double)N, sqrt(nrm2[0]/(double)N), nrmI[0],
1514: nrm1[1]/(double)N, sqrt(nrm2[1]/(double)N), nrmI[1],
1515: nrm1[2]/(double)N, sqrt(nrm2[2]/(double)N), nrmI[2] );
1516: */
1517: DMDestroy(&da_Stokes_analytic);
1518: VecDestroy(&X_analytic);
1519: }
1522: KSPDestroy(&ksp_S);
1523: VecDestroy(&X);
1524: VecDestroy(&f);
1525: MatDestroy(&A);
1526: MatDestroy(&B);
1528: DMDestroy(&da_Stokes);
1529: DMDestroy(&da_prop);
1531: VecDestroy(&properties);
1532: VecDestroy(&l_properties);
1534: return(0);
1535: }
1539: int main(int argc,char **args)
1540: {
1542: PetscInt mx,my;
1544: PetscInitialize(&argc,&args,(char *)0,help);
1546: mx = my = 10;
1547: PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,PETSC_NULL);
1548: PetscOptionsGetInt(PETSC_NULL,"-my",&my,PETSC_NULL);
1550: solve_stokes_2d_coupled(mx,my);
1552: PetscFinalize();
1553: return 0;
1554: }
1556: /* -------------------------- helpers for boundary conditions -------------------------------- */
1560: static PetscErrorCode BCApply_EAST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1561: {
1562: DM cda;
1563: Vec coords;
1564: PetscInt si,sj,nx,ny,i,j;
1565: PetscInt M,N;
1566: DMDACoor2d **_coords;
1567: PetscInt *g_idx;
1568: PetscInt *bc_global_ids;
1569: PetscScalar *bc_vals;
1570: PetscInt nbcs;
1571: PetscInt n_dofs;
1575: /* enforce bc's */
1576: DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);
1578: DMDAGetCoordinateDA(da,&cda);
1579: DMDAGetGhostedCoordinates(da,&coords);
1580: DMDAVecGetArray(cda,coords,&_coords);
1581: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1582: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1584: /* /// */
1586: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1587: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1589: /* init the entries to -1 so VecSetValues will ignore them */
1590: for (i = 0; i < ny*n_dofs; i++) {
1591: bc_global_ids[i] = -1;
1592: }
1594: i = nx-1;
1595: for (j = 0; j < ny; j++) {
1596: PetscInt local_id;
1598: local_id = i+j*nx;
1600: bc_global_ids[j] = g_idx[ n_dofs*local_id+d_idx ];
1602: bc_vals[j] = bc_val;
1603: }
1604: nbcs = 0;
1605: if ((si+nx) == (M)) {
1606: nbcs = ny;
1607: }
1609: if (b != PETSC_NULL) {
1610: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1611: VecAssemblyBegin(b);
1612: VecAssemblyEnd(b);
1613: }
1614: if (A != PETSC_NULL) {
1615: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1616: }
1619: PetscFree(bc_vals);
1620: PetscFree(bc_global_ids);
1622: DMDAVecRestoreArray(cda,coords,&_coords);
1623: return(0);
1624: }
1628: static PetscErrorCode BCApply_WEST(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1629: {
1630: DM cda;
1631: Vec coords;
1632: PetscInt si,sj,nx,ny,i,j;
1633: PetscInt M,N;
1634: DMDACoor2d **_coords;
1635: PetscInt *g_idx;
1636: PetscInt *bc_global_ids;
1637: PetscScalar *bc_vals;
1638: PetscInt nbcs;
1639: PetscInt n_dofs;
1643: /* enforce bc's */
1644: DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);
1646: DMDAGetCoordinateDA(da,&cda);
1647: DMDAGetGhostedCoordinates(da,&coords);
1648: DMDAVecGetArray(cda,coords,&_coords);
1649: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1650: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1652: /* /// */
1654: PetscMalloc(sizeof(PetscInt)*ny*n_dofs,&bc_global_ids);
1655: PetscMalloc(sizeof(PetscScalar)*ny*n_dofs,&bc_vals);
1657: /* init the entries to -1 so VecSetValues will ignore them */
1658: for (i = 0; i < ny*n_dofs; i++) {
1659: bc_global_ids[i] = -1;
1660: }
1662: i = 0;
1663: for (j = 0; j < ny; j++) {
1664: PetscInt local_id;
1666: local_id = i+j*nx;
1668: bc_global_ids[j] = g_idx[ n_dofs*local_id+d_idx ];
1670: bc_vals[j] = bc_val;
1671: }
1672: nbcs = 0;
1673: if (si == 0) {
1674: nbcs = ny;
1675: }
1677: if (b != PETSC_NULL) {
1678: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1679: VecAssemblyBegin(b);
1680: VecAssemblyEnd(b);
1681: }
1682: if (A != PETSC_NULL) {
1683: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1684: }
1687: PetscFree(bc_vals);
1688: PetscFree(bc_global_ids);
1690: DMDAVecRestoreArray(cda,coords,&_coords);
1691: return(0);
1692: }
1696: static PetscErrorCode BCApply_NORTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1697: {
1698: DM cda;
1699: Vec coords;
1700: PetscInt si,sj,nx,ny,i,j;
1701: PetscInt M,N;
1702: DMDACoor2d **_coords;
1703: PetscInt *g_idx;
1704: PetscInt *bc_global_ids;
1705: PetscScalar *bc_vals;
1706: PetscInt nbcs;
1707: PetscInt n_dofs;
1711: /* enforce bc's */
1712: DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);
1714: DMDAGetCoordinateDA(da,&cda);
1715: DMDAGetGhostedCoordinates(da,&coords);
1716: DMDAVecGetArray(cda,coords,&_coords);
1717: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1718: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1720: /* /// */
1722: PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1723: PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);
1725: /* init the entries to -1 so VecSetValues will ignore them */
1726: for (i = 0; i < nx; i++) {
1727: bc_global_ids[i] = -1;
1728: }
1730: j = ny-1;
1731: for (i = 0; i < nx; i++) {
1732: PetscInt local_id;
1734: local_id = i+j*nx;
1736: bc_global_ids[i] = g_idx[ n_dofs*local_id+d_idx ];
1738: bc_vals[i] = bc_val;
1739: }
1740: nbcs = 0;
1741: if ((sj+ny) == (N)) {
1742: nbcs = nx;
1743: }
1745: if (b != PETSC_NULL) {
1746: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1747: VecAssemblyBegin(b);
1748: VecAssemblyEnd(b);
1749: }
1750: if (A != PETSC_NULL) {
1751: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1752: }
1755: PetscFree(bc_vals);
1756: PetscFree(bc_global_ids);
1758: DMDAVecRestoreArray(cda,coords,&_coords);
1759: return(0);
1760: }
1764: static PetscErrorCode BCApply_SOUTH(DM da,PetscInt d_idx,PetscScalar bc_val,Mat A,Vec b)
1765: {
1766: DM cda;
1767: Vec coords;
1768: PetscInt si,sj,nx,ny,i,j;
1769: PetscInt M,N;
1770: DMDACoor2d **_coords;
1771: PetscInt *g_idx;
1772: PetscInt *bc_global_ids;
1773: PetscScalar *bc_vals;
1774: PetscInt nbcs;
1775: PetscInt n_dofs;
1779: /* enforce bc's */
1780: DMDAGetGlobalIndices(da,PETSC_NULL,&g_idx);
1782: DMDAGetCoordinateDA(da,&cda);
1783: DMDAGetGhostedCoordinates(da,&coords);
1784: DMDAVecGetArray(cda,coords,&_coords);
1785: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1786: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1788: /* /// */
1790: PetscMalloc(sizeof(PetscInt)*nx,&bc_global_ids);
1791: PetscMalloc(sizeof(PetscScalar)*nx,&bc_vals);
1793: /* init the entries to -1 so VecSetValues will ignore them */
1794: for (i = 0; i < nx; i++) {
1795: bc_global_ids[i] = -1;
1796: }
1798: j = 0;
1799: for (i = 0; i < nx; i++) {
1800: PetscInt local_id;
1802: local_id = i+j*nx;
1804: bc_global_ids[i] = g_idx[ n_dofs*local_id+d_idx ];
1806: bc_vals[i] = bc_val;
1807: }
1808: nbcs = 0;
1809: if (sj == 0) {
1810: nbcs = nx;
1811: }
1814: if (b != PETSC_NULL) {
1815: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1816: VecAssemblyBegin(b);
1817: VecAssemblyEnd(b);
1818: }
1819: if (A != PETSC_NULL) {
1820: MatZeroRows(A,nbcs,bc_global_ids,1.0,0,0);
1821: }
1824: PetscFree(bc_vals);
1825: PetscFree(bc_global_ids);
1827: DMDAVecRestoreArray(cda,coords,&_coords);
1828: return(0);
1829: }
1833: /*
1834: Free slip sides.
1835: */
1838: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1839: {
1843: BCApply_NORTH(da_Stokes,1,0.0,A,f);
1844: BCApply_EAST(da_Stokes,0,0.0,A,f);
1845: BCApply_SOUTH(da_Stokes,1,0.0,A,f);
1846: BCApply_WEST(da_Stokes,0,0.0,A,f);
1847: return(0);
1848: }