Actual source code: ex43.c
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: -glvis : Visualizes coefficients and u,p solution through GLVIs (use -viewer_glvis_dmda_bs 2,1 to visualize velocity as a vector)\n\
11: -c_str : Indicates the structure of the coefficients to use \n"
12: "\
13: -c_str 0 => Coefficient definition for an analytic solution with a vertical jump in viscosity at x = xc \n\
14: This problem is driven by the forcing function f(x,y) = (0, sin(nz pi y)cos(pi x) \n\
15: Parameters: \n\
16: -solcx_eta0 : Viscosity to the left of the interface \n\
17: -solcx_eta1 : Viscosity to the right of the interface \n\
18: -solcx_xc : Location of the interface \n\
19: -solcx_nz : Wavenumber in the y direction \n"
20: "\
21: -c_str 1 => Coefficient definition for a dense rectangular blob located at the center of the domain \n\
22: Parameters: \n\
23: -sinker_eta0 : Viscosity of the background fluid \n\
24: -sinker_eta1 : Viscosity of the blob \n\
25: -sinker_dx : Width of the blob \n\
26: -sinker_dy : Height of the blob \n"
27: "\
28: -c_str 2 => Coefficient definition for a dense circular blob located at the center of the domain \n\
29: Parameters: \n\
30: -sinker_eta0 : Viscosity of the background fluid \n\
31: -sinker_eta1 : Viscosity of the blob \n\
32: -sinker_r : Radius of the blob \n"
33: "\
34: -c_str 3 => Coefficient definition for a dense circular and rectangular inclusion (located at the center of the domain) \n\
35: -sinker_eta0 : Viscosity of the background fluid \n\
36: -sinker_eta1 : Viscosity of the two inclusions \n\
37: -sinker_r : Radius of the circular inclusion \n\
38: -sinker_c0x : Origin (x-coord) of the circular inclusion \n\
39: -sinker_c0y : Origin (y-coord) of the circular inclusion \n\
40: -sinker_dx : Width of the rectangular inclusion \n\
41: -sinker_dy : Height of the rectangular inclusion \n\
42: -sinker_phi : Rotation angle of the rectangular inclusion \n"
43: "\
44: -c_str 4 => Coefficient definition for checkerboard jumps aligned with the domain decomposition \n\
45: -jump_eta0 : Viscosity for black subdomains \n\
46: -jump_magnitude : Magnitude of jumps. White subdomains will have eta = eta0*10^magnitude \n\
47: -jump_nz : Wavenumber in the y direction for rhs \n"
48: "\
49: -use_gp_coords : Evaluate the viscosity and force term at the global coordinates of each quadrature point \n\
50: By default, the viscosity and force term are evaulated at the element center and applied as a constant over the entire element \n";
52: /* Contributed by Dave May */
54: #include <petscksp.h>
55: #include <petscdm.h>
56: #include <petscdmda.h>
58: /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
59: #include "ex43-solcx.h"
61: static PetscErrorCode DMDABCApplyFreeSlip(DM,Mat,Vec);
63: #define NSD 2 /* number of spatial dimensions */
64: #define NODES_PER_EL 4 /* nodes per element */
65: #define U_DOFS 2 /* degrees of freedom per velocity node */
66: #define P_DOFS 1 /* degrees of freedom per pressure node */
67: #define GAUSS_POINTS 4
69: /* Gauss point based evaluation 8+4+4+4 = 20 */
70: typedef struct {
71: PetscScalar gp_coords[2*GAUSS_POINTS];
72: PetscScalar eta[GAUSS_POINTS];
73: PetscScalar fx[GAUSS_POINTS];
74: PetscScalar fy[GAUSS_POINTS];
75: } GaussPointCoefficients;
77: typedef struct {
78: PetscScalar u_dof;
79: PetscScalar v_dof;
80: PetscScalar p_dof;
81: } StokesDOF;
83: static PetscErrorCode glvis_extract_eta(PetscObject oV,PetscInt nf, PetscObject oVf[], void *ctx)
84: {
85: DM properties_da = (DM)(ctx),stokes_da;
86: Vec V = (Vec)oV, *Vf = (Vec*)oVf;
87: GaussPointCoefficients **props;
88: PetscInt sex,sey,mx,my;
89: PetscInt ei,ej,p,cum;
90: PetscScalar *array;
93: VecGetDM(Vf[0],&stokes_da);
94: DMDAVecGetArrayRead(properties_da,V,&props);
95: DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
96: DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
97: VecGetArray(Vf[0],&array);
98: cum = 0;
99: for (ej = sey; ej < sey+my; ej++) {
100: for (ei = sex; ei < sex+mx; ei++) {
101: for (p = 0; p < GAUSS_POINTS; p++) {
102: array[cum++] = props[ej][ei].eta[p];
103: }
104: }
105: }
106: VecRestoreArray(Vf[0],&array);
107: DMDAVecRestoreArrayRead(properties_da,V,&props);
108: return 0;
109: }
111: /*
112: Element: Local basis function ordering
113: 1-----2
114: | |
115: | |
116: 0-----3
117: */
118: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
119: {
120: PetscScalar xi = _xi[0];
121: PetscScalar eta = _xi[1];
123: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
124: Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
125: Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
126: Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
127: }
129: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
130: {
131: PetscScalar xi = _xi[0];
132: PetscScalar eta = _xi[1];
134: GNi[0][0] = -0.25*(1.0-eta);
135: GNi[0][1] = -0.25*(1.0+eta);
136: GNi[0][2] = 0.25*(1.0+eta);
137: GNi[0][3] = 0.25*(1.0-eta);
139: GNi[1][0] = -0.25*(1.0-xi);
140: GNi[1][1] = 0.25*(1.0-xi);
141: GNi[1][2] = 0.25*(1.0+xi);
142: GNi[1][3] = -0.25*(1.0+xi);
143: }
145: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
146: {
147: PetscScalar J00,J01,J10,J11,J;
148: PetscScalar iJ00,iJ01,iJ10,iJ11;
149: PetscInt i;
151: J00 = J01 = J10 = J11 = 0.0;
152: for (i = 0; i < NODES_PER_EL; i++) {
153: PetscScalar cx = coords[2*i];
154: PetscScalar cy = coords[2*i+1];
156: J00 = J00+GNi[0][i]*cx; /* J_xx = dx/dxi */
157: J01 = J01+GNi[0][i]*cy; /* J_xy = dy/dxi */
158: J10 = J10+GNi[1][i]*cx; /* J_yx = dx/deta */
159: J11 = J11+GNi[1][i]*cy; /* J_yy = dy/deta */
160: }
161: J = (J00*J11)-(J01*J10);
163: iJ00 = J11/J;
164: iJ01 = -J01/J;
165: iJ10 = -J10/J;
166: iJ11 = J00/J;
168: for (i = 0; i < NODES_PER_EL; i++) {
169: GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
170: GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
171: }
173: if (det_J) *det_J = J;
174: }
176: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
177: {
178: *ngp = 4;
179: gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919;
180: gp_xi[1][0] = -0.57735026919; gp_xi[1][1] = 0.57735026919;
181: gp_xi[2][0] = 0.57735026919; gp_xi[2][1] = 0.57735026919;
182: gp_xi[3][0] = 0.57735026919; gp_xi[3][1] = -0.57735026919;
183: gp_weight[0] = 1.0;
184: gp_weight[1] = 1.0;
185: gp_weight[2] = 1.0;
186: gp_weight[3] = 1.0;
187: }
189: /*
190: i,j are the element indices
191: The unknown is a vector quantity.
192: The s[].c is used to indicate the degree of freedom.
193: */
194: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
195: {
197: /* velocity */
198: /* node 0 */
199: s_u[0].i = i; s_u[0].j = j; s_u[0].c = 0; /* Vx0 */
200: s_u[1].i = i; s_u[1].j = j; s_u[1].c = 1; /* Vy0 */
202: /* node 1 */
203: s_u[2].i = i; s_u[2].j = j+1; s_u[2].c = 0; /* Vx1 */
204: s_u[3].i = i; s_u[3].j = j+1; s_u[3].c = 1; /* Vy1 */
206: /* node 2 */
207: s_u[4].i = i+1; s_u[4].j = j+1; s_u[4].c = 0; /* Vx2 */
208: s_u[5].i = i+1; s_u[5].j = j+1; s_u[5].c = 1; /* Vy2 */
210: /* node 3 */
211: s_u[6].i = i+1; s_u[6].j = j; s_u[6].c = 0; /* Vx3 */
212: s_u[7].i = i+1; s_u[7].j = j; s_u[7].c = 1; /* Vy3 */
214: /* pressure */
215: s_p[0].i = i; s_p[0].j = j; s_p[0].c = 2; /* P0 */
216: s_p[1].i = i; s_p[1].j = j+1; s_p[1].c = 2; /* P1 */
217: s_p[2].i = i+1; s_p[2].j = j+1; s_p[2].c = 2; /* P2 */
218: s_p[3].i = i+1; s_p[3].j = j; s_p[3].c = 2; /* P3 */
219: return 0;
220: }
222: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
223: {
224: PetscMPIInt rank;
225: PetscInt proc_I,proc_J;
226: PetscInt cpu_x,cpu_y;
227: PetscInt local_mx,local_my;
228: Vec vlx,vly;
229: PetscInt *LX,*LY,i;
230: PetscScalar *_a;
231: Vec V_SEQ;
232: VecScatter ctx;
235: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
237: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
239: proc_J = rank/cpu_x;
240: proc_I = rank-cpu_x*proc_J;
242: PetscMalloc1(cpu_x,&LX);
243: PetscMalloc1(cpu_y,&LY);
245: DMDAGetElementsSizes(da,&local_mx,&local_my,NULL);
246: VecCreate(PETSC_COMM_WORLD,&vlx);
247: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
248: VecSetFromOptions(vlx);
250: VecCreate(PETSC_COMM_WORLD,&vly);
251: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
252: VecSetFromOptions(vly);
254: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
255: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
256: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
257: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
259: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
260: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
261: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
262: VecGetArray(V_SEQ,&_a);
263: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
264: VecRestoreArray(V_SEQ,&_a);
265: VecScatterDestroy(&ctx);
266: VecDestroy(&V_SEQ);
268: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
269: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
270: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
271: VecGetArray(V_SEQ,&_a);
272: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
273: VecRestoreArray(V_SEQ,&_a);
274: VecScatterDestroy(&ctx);
275: VecDestroy(&V_SEQ);
277: *_lx = LX;
278: *_ly = LY;
280: VecDestroy(&vlx);
281: VecDestroy(&vly);
282: return 0;
283: }
285: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
286: {
287: DM cda;
288: Vec coords;
289: DMDACoor2d **_coords;
290: PetscInt si,sj,nx,ny,i,j;
291: FILE *fp;
292: char fname[PETSC_MAX_PATH_LEN];
293: PetscMPIInt rank;
296: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
297: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
298: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
301: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
303: DMGetCoordinateDM(da,&cda);
304: DMGetCoordinatesLocal(da,&coords);
305: DMDAVecGetArrayRead(cda,coords,&_coords);
306: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
307: for (j = sj; j < sj+ny-1; j++) {
308: for (i = si; i < si+nx-1; i++) {
309: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
310: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i].x),(double)PetscRealPart(_coords[j+1][i].y));
311: 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));
312: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i+1].x),(double)PetscRealPart(_coords[j][i+1].y));
313: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
314: }
315: }
316: DMDAVecRestoreArrayRead(cda,coords,&_coords);
318: PetscFClose(PETSC_COMM_SELF,fp);
319: return 0;
320: }
322: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
323: {
324: DM cda;
325: Vec coords,local_fields;
326: DMDACoor2d **_coords;
327: FILE *fp;
328: char fname[PETSC_MAX_PATH_LEN];
329: PetscMPIInt rank;
330: PetscInt si,sj,nx,ny,i,j;
331: PetscInt n_dofs,d;
332: PetscScalar *_fields;
335: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
336: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
337: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
340: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
341: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
342: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
343: for (d = 0; d < n_dofs; d++) {
344: const char *field_name;
345: DMDAGetFieldName(da,d,&field_name);
346: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
347: }
348: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
350: DMGetCoordinateDM(da,&cda);
351: DMGetCoordinatesLocal(da,&coords);
352: DMDAVecGetArray(cda,coords,&_coords);
353: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
355: DMCreateLocalVector(da,&local_fields);
356: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
357: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
358: VecGetArray(local_fields,&_fields);
360: for (j = sj; j < sj+ny; j++) {
361: for (i = si; i < si+nx; i++) {
362: PetscScalar coord_x,coord_y;
363: PetscScalar field_d;
365: coord_x = _coords[j][i].x;
366: coord_y = _coords[j][i].y;
368: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
369: for (d = 0; d < n_dofs; d++) {
370: field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
371: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
372: }
373: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
374: }
375: }
376: VecRestoreArray(local_fields,&_fields);
377: VecDestroy(&local_fields);
379: DMDAVecRestoreArray(cda,coords,&_coords);
381: PetscFClose(PETSC_COMM_SELF,fp);
382: return 0;
383: }
385: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
386: {
387: DM cda;
388: Vec local_fields;
389: FILE *fp;
390: char fname[PETSC_MAX_PATH_LEN];
391: PetscMPIInt rank;
392: PetscInt si,sj,nx,ny,i,j,p;
393: PetscInt n_dofs,d;
394: GaussPointCoefficients **_coefficients;
397: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
398: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
399: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
402: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
403: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
404: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
405: for (d = 0; d < n_dofs; d++) {
406: const char *field_name;
407: DMDAGetFieldName(da,d,&field_name);
408: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
409: }
410: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
412: DMGetCoordinateDM(da,&cda);
413: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
415: DMCreateLocalVector(da,&local_fields);
416: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
417: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
418: DMDAVecGetArray(da,local_fields,&_coefficients);
420: for (j = sj; j < sj+ny; j++) {
421: for (i = si; i < si+nx; i++) {
422: PetscScalar coord_x,coord_y;
424: for (p = 0; p < GAUSS_POINTS; p++) {
425: coord_x = _coefficients[j][i].gp_coords[2*p];
426: coord_y = _coefficients[j][i].gp_coords[2*p+1];
428: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
430: 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]));
431: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
432: }
433: }
434: }
435: DMDAVecRestoreArray(da,local_fields,&_coefficients);
436: VecDestroy(&local_fields);
438: PetscFClose(PETSC_COMM_SELF,fp);
439: return 0;
440: }
442: 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)
443: {
444: PetscInt ij;
445: PetscInt r,c,nc;
447: nc = u_NPE*u_dof;
448: r = w_dof*wi+wd;
449: c = u_dof*ui+ud;
450: ij = r*nc+c;
451: return ij;
452: }
454: /*
455: D = [ 2.eta 0 0 ]
456: [ 0 2.eta 0 ]
457: [ 0 0 eta ]
459: B = [ d_dx 0 ]
460: [ 0 d_dy ]
461: [ d_dy d_dx ]
462: */
463: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
464: {
465: PetscInt ngp;
466: PetscScalar gp_xi[GAUSS_POINTS][2];
467: PetscScalar gp_weight[GAUSS_POINTS];
468: PetscInt p,i,j,k;
469: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
470: PetscScalar J_p,tildeD[3];
471: PetscScalar B[3][U_DOFS*NODES_PER_EL];
473: /* define quadrature rule */
474: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
476: /* evaluate integral */
477: for (p = 0; p < ngp; p++) {
478: ConstructQ12D_GNi(gp_xi[p],GNi_p);
479: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
481: for (i = 0; i < NODES_PER_EL; i++) {
482: PetscScalar d_dx_i = GNx_p[0][i];
483: PetscScalar d_dy_i = GNx_p[1][i];
485: B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
486: B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
487: B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
488: }
490: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
491: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
492: tildeD[2] = gp_weight[p]*J_p*eta[p];
494: /* form Bt tildeD B */
495: /*
496: Ke_ij = Bt_ik . D_kl . B_lj
497: = B_ki . D_kl . B_lj
498: = B_ki . D_kk . B_kj
499: */
500: for (i = 0; i < 8; i++) {
501: for (j = 0; j < 8; j++) {
502: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
503: Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
504: }
505: }
506: }
507: }
508: }
510: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
511: {
512: PetscInt ngp;
513: PetscScalar gp_xi[GAUSS_POINTS][2];
514: PetscScalar gp_weight[GAUSS_POINTS];
515: PetscInt p,i,j,di;
516: PetscScalar Ni_p[NODES_PER_EL];
517: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
518: PetscScalar J_p,fac;
520: /* define quadrature rule */
521: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
523: /* evaluate integral */
524: for (p = 0; p < ngp; p++) {
525: ConstructQ12D_Ni(gp_xi[p],Ni_p);
526: ConstructQ12D_GNi(gp_xi[p],GNi_p);
527: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
528: fac = gp_weight[p]*J_p;
530: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
531: for (di = 0; di < NSD; di++) { /* u dofs */
532: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
533: PetscInt IJ;
534: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);
536: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
537: }
538: }
539: }
540: }
541: }
543: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
544: {
545: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
546: PetscInt i,j;
547: PetscInt nr_g,nc_g;
549: PetscMemzero(Ge,sizeof(Ge));
550: FormGradientOperatorQ1(Ge,coords);
552: nr_g = U_DOFS*NODES_PER_EL;
553: nc_g = P_DOFS*NODES_PER_EL;
555: for (i = 0; i < nr_g; i++) {
556: for (j = 0; j < nc_g; j++) {
557: De[nr_g*j+i] = Ge[nc_g*i+j];
558: }
559: }
560: }
562: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
563: {
564: PetscInt ngp;
565: PetscScalar gp_xi[GAUSS_POINTS][2];
566: PetscScalar gp_weight[GAUSS_POINTS];
567: PetscInt p,i,j;
568: PetscScalar Ni_p[NODES_PER_EL];
569: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
570: PetscScalar J_p,fac,eta_avg;
572: /* define quadrature rule */
573: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
575: /* evaluate integral */
576: for (p = 0; p < ngp; p++) {
577: ConstructQ12D_Ni(gp_xi[p],Ni_p);
578: ConstructQ12D_GNi(gp_xi[p],GNi_p);
579: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
580: fac = gp_weight[p]*J_p;
582: for (i = 0; i < NODES_PER_EL; i++) {
583: for (j = 0; j < NODES_PER_EL; j++) {
584: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
585: }
586: }
587: }
589: /* scale */
590: eta_avg = 0.0;
591: for (p = 0; p < ngp; p++) eta_avg += eta[p];
592: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
593: fac = 1.0/eta_avg;
594: for (i = 0; i < NODES_PER_EL; i++) {
595: for (j = 0; j < NODES_PER_EL; j++) {
596: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
597: }
598: }
599: }
601: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
602: {
603: PetscInt ngp;
604: PetscScalar gp_xi[GAUSS_POINTS][2];
605: PetscScalar gp_weight[GAUSS_POINTS];
606: PetscInt p,i,j;
607: PetscScalar Ni_p[NODES_PER_EL];
608: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
609: PetscScalar J_p,fac,eta_avg;
611: /* define quadrature rule */
612: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
614: /* evaluate integral */
615: for (p = 0; p < ngp; p++) {
616: ConstructQ12D_Ni(gp_xi[p],Ni_p);
617: ConstructQ12D_GNi(gp_xi[p],GNi_p);
618: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
619: fac = gp_weight[p]*J_p;
621: for (i = 0; i < NODES_PER_EL; i++) {
622: for (j = 0; j < NODES_PER_EL; j++) {
623: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
624: }
625: }
626: }
628: /* scale */
629: eta_avg = 0.0;
630: for (p = 0; p < ngp; p++) eta_avg += eta[p];
631: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
632: fac = 1.0/eta_avg;
633: for (i = 0; i < NODES_PER_EL; i++) {
634: for (j = 0; j < NODES_PER_EL; j++) {
635: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
636: }
637: }
638: }
640: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
641: {
642: PetscInt ngp;
643: PetscScalar gp_xi[GAUSS_POINTS][2];
644: PetscScalar gp_weight[GAUSS_POINTS];
645: PetscInt p,i;
646: PetscScalar Ni_p[NODES_PER_EL];
647: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
648: PetscScalar J_p,fac;
650: /* define quadrature rule */
651: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
653: /* evaluate integral */
654: for (p = 0; p < ngp; p++) {
655: ConstructQ12D_Ni(gp_xi[p],Ni_p);
656: ConstructQ12D_GNi(gp_xi[p],GNi_p);
657: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
658: fac = gp_weight[p]*J_p;
660: for (i = 0; i < NODES_PER_EL; i++) {
661: Fe[NSD*i] += fac*Ni_p[i]*fx[p];
662: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
663: }
664: }
665: }
667: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
668: {
670: /* get coords for the element */
671: el_coords[NSD*0] = _coords[ej][ei].x; el_coords[NSD*0+1] = _coords[ej][ei].y;
672: el_coords[NSD*1] = _coords[ej+1][ei].x; el_coords[NSD*1+1] = _coords[ej+1][ei].y;
673: el_coords[NSD*2] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
674: el_coords[NSD*3] = _coords[ej][ei+1].x; el_coords[NSD*3+1] = _coords[ej][ei+1].y;
675: return 0;
676: }
678: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
679: {
680: DM cda;
681: Vec coords;
682: DMDACoor2d **_coords;
683: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
684: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
685: PetscInt sex,sey,mx,my;
686: PetscInt ei,ej;
687: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
688: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
689: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
690: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
691: PetscScalar el_coords[NODES_PER_EL*NSD];
692: Vec local_properties;
693: GaussPointCoefficients **props;
694: PetscScalar *prop_eta;
697: /* setup for coords */
698: DMGetCoordinateDM(stokes_da,&cda);
699: DMGetCoordinatesLocal(stokes_da,&coords);
700: DMDAVecGetArray(cda,coords,&_coords);
702: /* setup for coefficients */
703: DMCreateLocalVector(properties_da,&local_properties);
704: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
705: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
706: DMDAVecGetArray(properties_da,local_properties,&props);
708: DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
709: DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
710: for (ej = sey; ej < sey+my; ej++) {
711: for (ei = sex; ei < sex+mx; ei++) {
712: /* get coords for the element */
713: GetElementCoords(_coords,ei,ej,el_coords);
715: /* get coefficients for the element */
716: prop_eta = props[ej][ei].eta;
718: /* initialise element stiffness matrix */
719: PetscMemzero(Ae,sizeof(Ae));
720: PetscMemzero(Ge,sizeof(Ge));
721: PetscMemzero(De,sizeof(De));
722: PetscMemzero(Ce,sizeof(Ce));
724: /* form element stiffness matrix */
725: FormStressOperatorQ1(Ae,el_coords,prop_eta);
726: FormGradientOperatorQ1(Ge,el_coords);
727: FormDivergenceOperatorQ1(De,el_coords);
728: FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);
730: /* insert element matrix into global matrix */
731: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
732: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
733: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
734: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
735: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
736: }
737: }
738: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
739: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
741: DMDAVecRestoreArray(cda,coords,&_coords);
743: DMDAVecRestoreArray(properties_da,local_properties,&props);
744: VecDestroy(&local_properties);
745: return 0;
746: }
748: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
749: {
750: DM cda;
751: Vec coords;
752: DMDACoor2d **_coords;
753: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
754: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
755: PetscInt sex,sey,mx,my;
756: PetscInt ei,ej;
757: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
758: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
759: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
760: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
761: PetscScalar el_coords[NODES_PER_EL*NSD];
762: Vec local_properties;
763: GaussPointCoefficients **props;
764: PetscScalar *prop_eta;
767: /* setup for coords */
768: DMGetCoordinateDM(stokes_da,&cda);
769: DMGetCoordinatesLocal(stokes_da,&coords);
770: DMDAVecGetArray(cda,coords,&_coords);
772: /* setup for coefficients */
773: DMCreateLocalVector(properties_da,&local_properties);
774: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
775: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
776: DMDAVecGetArray(properties_da,local_properties,&props);
778: DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
779: DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
780: for (ej = sey; ej < sey+my; ej++) {
781: for (ei = sex; ei < sex+mx; ei++) {
782: /* get coords for the element */
783: GetElementCoords(_coords,ei,ej,el_coords);
785: /* get coefficients for the element */
786: prop_eta = props[ej][ei].eta;
788: /* initialise element stiffness matrix */
789: PetscMemzero(Ae,sizeof(Ae));
790: PetscMemzero(Ge,sizeof(Ge));
791: PetscMemzero(De,sizeof(De));
792: PetscMemzero(Ce,sizeof(Ce));
794: /* form element stiffness matrix */
795: FormStressOperatorQ1(Ae,el_coords,prop_eta);
796: FormGradientOperatorQ1(Ge,el_coords);
797: FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);
799: /* insert element matrix into global matrix */
800: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
801: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
802: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
803: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
804: }
805: }
806: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
807: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
809: DMDAVecRestoreArray(cda,coords,&_coords);
811: DMDAVecRestoreArray(properties_da,local_properties,&props);
812: VecDestroy(&local_properties);
813: return 0;
814: }
816: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
817: {
818: PetscInt n;
821: for (n = 0; n < 4; n++) {
822: 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];
823: 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];
824: 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];
825: }
826: return 0;
827: }
829: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
830: {
831: DM cda;
832: Vec coords;
833: DMDACoor2d **_coords;
834: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
835: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
836: PetscInt sex,sey,mx,my;
837: PetscInt ei,ej;
838: PetscScalar Fe[NODES_PER_EL*U_DOFS];
839: PetscScalar He[NODES_PER_EL*P_DOFS];
840: PetscScalar el_coords[NODES_PER_EL*NSD];
841: Vec local_properties;
842: GaussPointCoefficients **props;
843: PetscScalar *prop_fx,*prop_fy;
844: Vec local_F;
845: StokesDOF **ff;
848: /* setup for coords */
849: DMGetCoordinateDM(stokes_da,&cda);
850: DMGetCoordinatesLocal(stokes_da,&coords);
851: DMDAVecGetArray(cda,coords,&_coords);
853: /* setup for coefficients */
854: DMGetLocalVector(properties_da,&local_properties);
855: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
856: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
857: DMDAVecGetArray(properties_da,local_properties,&props);
859: /* get access to the vector */
860: DMGetLocalVector(stokes_da,&local_F);
861: VecZeroEntries(local_F);
862: DMDAVecGetArray(stokes_da,local_F,&ff);
864: DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
865: DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
866: for (ej = sey; ej < sey+my; ej++) {
867: for (ei = sex; ei < sex+mx; ei++) {
868: /* get coords for the element */
869: GetElementCoords(_coords,ei,ej,el_coords);
871: /* get coefficients for the element */
872: prop_fx = props[ej][ei].fx;
873: prop_fy = props[ej][ei].fy;
875: /* initialise element stiffness matrix */
876: PetscMemzero(Fe,sizeof(Fe));
877: PetscMemzero(He,sizeof(He));
879: /* form element stiffness matrix */
880: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
882: /* insert element matrix into global matrix */
883: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
885: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
886: }
887: }
889: DMDAVecRestoreArray(stokes_da,local_F,&ff);
890: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
891: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
892: DMRestoreLocalVector(stokes_da,&local_F);
894: DMDAVecRestoreArray(cda,coords,&_coords);
896: DMDAVecRestoreArray(properties_da,local_properties,&props);
897: DMRestoreLocalVector(properties_da,&local_properties);
898: return 0;
899: }
901: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM *_da,Vec *_X)
902: {
903: DM da,cda;
904: Vec X;
905: StokesDOF **_stokes;
906: Vec coords;
907: DMDACoor2d **_coords;
908: PetscInt si,sj,ei,ej,i,j;
911: 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);
912: DMSetFromOptions(da);
913: DMSetUp(da);
914: DMDASetFieldName(da,0,"anlytic_Vx");
915: DMDASetFieldName(da,1,"anlytic_Vy");
916: DMDASetFieldName(da,2,"analytic_P");
918: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.,0.);
920: DMGetCoordinatesLocal(da,&coords);
921: DMGetCoordinateDM(da,&cda);
922: DMDAVecGetArray(cda,coords,&_coords);
924: DMCreateGlobalVector(da,&X);
925: DMDAVecGetArray(da,X,&_stokes);
927: DMDAGetCorners(da,&si,&sj,0,&ei,&ej,0);
928: for (j = sj; j < sj+ej; j++) {
929: for (i = si; i < si+ei; i++) {
930: PetscReal pos[2],pressure,vel[2],total_stress[3],strain_rate[3];
932: pos[0] = PetscRealPart(_coords[j][i].x);
933: pos[1] = PetscRealPart(_coords[j][i].y);
935: evaluate_solCx(pos,eta0,eta1,xc,nz,vel,&pressure,total_stress,strain_rate);
937: _stokes[j][i].u_dof = vel[0];
938: _stokes[j][i].v_dof = vel[1];
939: _stokes[j][i].p_dof = pressure;
940: }
941: }
942: DMDAVecRestoreArray(da,X,&_stokes);
943: DMDAVecRestoreArray(cda,coords,&_coords);
945: *_da = da;
946: *_X = X;
947: return 0;
948: }
950: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
951: {
953: /* get the nodal fields */
954: 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;
955: 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;
956: 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;
957: 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;
958: return 0;
959: }
961: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
962: {
963: DM cda;
964: Vec coords,X_analytic_local,X_local;
965: DMDACoor2d **_coords;
966: PetscInt sex,sey,mx,my;
967: PetscInt ei,ej;
968: PetscScalar el_coords[NODES_PER_EL*NSD];
969: StokesDOF **stokes_analytic,**stokes;
970: StokesDOF stokes_analytic_e[4],stokes_e[4];
972: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
973: PetscScalar Ni_p[NODES_PER_EL];
974: PetscInt ngp;
975: PetscScalar gp_xi[GAUSS_POINTS][2];
976: PetscScalar gp_weight[GAUSS_POINTS];
977: PetscInt p,i;
978: PetscScalar J_p,fac;
979: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
980: PetscInt M;
981: PetscReal xymin[2],xymax[2];
984: /* define quadrature rule */
985: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
987: /* setup for coords */
988: DMGetCoordinateDM(stokes_da,&cda);
989: DMGetCoordinatesLocal(stokes_da,&coords);
990: DMDAVecGetArray(cda,coords,&_coords);
992: /* setup for analytic */
993: DMCreateLocalVector(stokes_da,&X_analytic_local);
994: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
995: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
996: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
998: /* setup for solution */
999: DMCreateLocalVector(stokes_da,&X_local);
1000: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1001: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1002: DMDAVecGetArray(stokes_da,X_local,&stokes);
1004: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1005: DMGetBoundingBox(stokes_da,xymin,xymax);
1007: h = (xymax[0]-xymin[0])/((PetscReal)M);
1009: tp_L2 = tu_L2 = tu_H1 = 0.0;
1011: DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
1012: DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
1013: for (ej = sey; ej < sey+my; ej++) {
1014: for (ei = sex; ei < sex+mx; ei++) {
1015: /* get coords for the element */
1016: GetElementCoords(_coords,ei,ej,el_coords);
1017: StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1018: StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);
1020: /* evaluate integral */
1021: p_e_L2 = 0.0;
1022: u_e_L2 = 0.0;
1023: u_e_H1 = 0.0;
1024: for (p = 0; p < ngp; p++) {
1025: ConstructQ12D_Ni(gp_xi[p],Ni_p);
1026: ConstructQ12D_GNi(gp_xi[p],GNi_p);
1027: ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1028: fac = gp_weight[p]*J_p;
1030: for (i = 0; i < NODES_PER_EL; i++) {
1031: PetscScalar u_error,v_error;
1033: 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);
1035: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1036: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1037: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);
1039: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1040: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error /* du/dy */
1041: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1042: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error); /* dv/dy */
1043: }
1044: }
1046: tp_L2 += p_e_L2;
1047: tu_L2 += u_e_L2;
1048: tu_H1 += u_e_H1;
1049: }
1050: }
1051: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1052: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1053: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1054: p_L2 = PetscSqrtScalar(p_L2);
1055: u_L2 = PetscSqrtScalar(u_L2);
1056: u_H1 = PetscSqrtScalar(u_H1);
1058: 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));
1060: DMDAVecRestoreArray(cda,coords,&_coords);
1062: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1063: VecDestroy(&X_analytic_local);
1064: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1065: VecDestroy(&X_local);
1066: return 0;
1067: }
1069: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1070: {
1071: DM da_Stokes,da_prop;
1072: PetscInt u_dof,p_dof,dof,stencil_width;
1073: Mat A,B;
1074: DM prop_cda,vel_cda;
1075: Vec prop_coords,vel_coords;
1076: PetscInt si,sj,nx,ny,i,j,p;
1077: Vec f,X;
1078: PetscInt prop_dof,prop_stencil_width;
1079: Vec properties,l_properties;
1080: PetscReal dx,dy;
1081: PetscInt M,N;
1082: DMDACoor2d **_prop_coords,**_vel_coords;
1083: GaussPointCoefficients **element_props;
1084: PetscInt its;
1085: KSP ksp_S;
1086: PetscInt coefficient_structure = 0;
1087: PetscInt cpu_x,cpu_y,*lx = NULL,*ly = NULL;
1088: PetscBool use_gp_coords = PETSC_FALSE,set,output_gnuplot = PETSC_FALSE,glvis = PETSC_FALSE,change = PETSC_FALSE;
1089: char filename[PETSC_MAX_PATH_LEN];
1093: PetscOptionsGetBool(NULL,NULL,"-gnuplot",&output_gnuplot,NULL);
1094: PetscOptionsGetBool(NULL,NULL,"-glvis",&glvis,NULL);
1095: PetscOptionsGetBool(NULL,NULL,"-change",&change,NULL);
1097: /* Generate the da for velocity and pressure */
1098: /*
1099: We use Q1 elements for the temperature.
1100: FEM has a 9-point stencil (BOX) or connectivity pattern
1101: Num nodes in each direction is mx+1, my+1
1102: */
1103: u_dof = U_DOFS; /* Vx, Vy - velocities */
1104: p_dof = P_DOFS; /* p - pressure */
1105: dof = u_dof+p_dof;
1106: stencil_width = 1;
1107: 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);
1109: DMSetMatType(da_Stokes,MATAIJ);
1110: DMSetFromOptions(da_Stokes);
1111: DMSetUp(da_Stokes);
1112: DMDASetFieldName(da_Stokes,0,"Vx");
1113: DMDASetFieldName(da_Stokes,1,"Vy");
1114: DMDASetFieldName(da_Stokes,2,"P");
1116: /* unit box [0,1] x [0,1] */
1117: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.,0.);
1119: /* Generate element properties, we will assume all material properties are constant over the element */
1120: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1121: DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1122: DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);
1124: prop_dof = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1125: prop_stencil_width = 0;
1126: 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);
1127: DMSetFromOptions(da_prop);
1128: DMSetUp(da_prop);
1129: PetscFree(lx);
1130: PetscFree(ly);
1132: /* define centroid positions */
1133: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1134: dx = 1.0/((PetscReal)(M));
1135: dy = 1.0/((PetscReal)(N));
1137: 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);
1139: /* define coefficients */
1140: PetscOptionsGetInt(NULL,NULL,"-c_str",&coefficient_structure,NULL);
1142: DMCreateGlobalVector(da_prop,&properties);
1143: DMCreateLocalVector(da_prop,&l_properties);
1144: DMDAVecGetArray(da_prop,l_properties,&element_props);
1146: DMGetCoordinateDM(da_prop,&prop_cda);
1147: DMGetCoordinatesLocal(da_prop,&prop_coords);
1148: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
1150: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
1152: DMGetCoordinateDM(da_Stokes,&vel_cda);
1153: DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1154: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1156: /* interpolate the coordinates */
1157: for (j = sj; j < sj+ny; j++) {
1158: for (i = si; i < si+nx; i++) {
1159: PetscInt ngp;
1160: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1161: PetscScalar el_coords[8];
1163: GetElementCoords(_vel_coords,i,j,el_coords);
1164: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1166: for (p = 0; p < GAUSS_POINTS; p++) {
1167: PetscScalar gp_x,gp_y;
1168: PetscInt n;
1169: PetscScalar xi_p[2],Ni_p[4];
1171: xi_p[0] = gp_xi[p][0];
1172: xi_p[1] = gp_xi[p][1];
1173: ConstructQ12D_Ni(xi_p,Ni_p);
1175: gp_x = 0.0;
1176: gp_y = 0.0;
1177: for (n = 0; n < NODES_PER_EL; n++) {
1178: gp_x = gp_x+Ni_p[n]*el_coords[2*n];
1179: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1180: }
1181: element_props[j][i].gp_coords[2*p] = gp_x;
1182: element_props[j][i].gp_coords[2*p+1] = gp_y;
1183: }
1184: }
1185: }
1187: /* define the coefficients */
1188: PetscOptionsGetBool(NULL,NULL,"-use_gp_coords",&use_gp_coords,NULL);
1190: for (j = sj; j < sj+ny; j++) {
1191: for (i = si; i < si+nx; i++) {
1192: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1193: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1194: PetscReal coord_x,coord_y;
1196: if (coefficient_structure == 0) {
1197: PetscReal opts_eta0,opts_eta1,opts_xc;
1198: PetscInt opts_nz;
1200: opts_eta0 = 1.0;
1201: opts_eta1 = 1.0;
1202: opts_xc = 0.5;
1203: opts_nz = 1;
1205: PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1206: PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1207: PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1208: PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);
1210: for (p = 0; p < GAUSS_POINTS; p++) {
1211: coord_x = centroid_x;
1212: coord_y = centroid_y;
1213: if (use_gp_coords) {
1214: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1215: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1216: }
1218: element_props[j][i].eta[p] = opts_eta0;
1219: if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;
1221: element_props[j][i].fx[p] = 0.0;
1222: element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1223: }
1224: } else if (coefficient_structure == 1) { /* square sinker */
1225: PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;
1227: opts_eta0 = 1.0;
1228: opts_eta1 = 1.0;
1229: opts_dx = 0.50;
1230: opts_dy = 0.50;
1232: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1233: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1234: PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1235: PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
1237: for (p = 0; p < GAUSS_POINTS; p++) {
1238: coord_x = centroid_x;
1239: coord_y = centroid_y;
1240: if (use_gp_coords) {
1241: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1242: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1243: }
1245: element_props[j][i].eta[p] = opts_eta0;
1246: element_props[j][i].fx[p] = 0.0;
1247: element_props[j][i].fy[p] = 0.0;
1249: if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1250: if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1251: element_props[j][i].eta[p] = opts_eta1;
1252: element_props[j][i].fx[p] = 0.0;
1253: element_props[j][i].fy[p] = -1.0;
1254: }
1255: }
1256: }
1257: } else if (coefficient_structure == 2) { /* circular sinker */
1258: PetscReal opts_eta0,opts_eta1,opts_r,radius2;
1260: opts_eta0 = 1.0;
1261: opts_eta1 = 1.0;
1262: opts_r = 0.25;
1264: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1265: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1266: PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);
1268: for (p = 0; p < GAUSS_POINTS; p++) {
1269: coord_x = centroid_x;
1270: coord_y = centroid_y;
1271: if (use_gp_coords) {
1272: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1273: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1274: }
1276: element_props[j][i].eta[p] = opts_eta0;
1277: element_props[j][i].fx[p] = 0.0;
1278: element_props[j][i].fy[p] = 0.0;
1280: radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1281: if (radius2 < opts_r*opts_r) {
1282: element_props[j][i].eta[p] = opts_eta1;
1283: element_props[j][i].fx[p] = 0.0;
1284: element_props[j][i].fy[p] = -1.0;
1285: }
1286: }
1287: } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1288: PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;
1290: opts_eta0 = 1.0;
1291: opts_eta1 = 1.0;
1292: opts_r = 0.25;
1293: opts_c0x = 0.35; /* circle center */
1294: opts_c0y = 0.35;
1295: opts_s0x = 0.7; /* square center */
1296: opts_s0y = 0.7;
1297: opts_dx = 0.25;
1298: opts_dy = 0.25;
1299: opts_phi = 25;
1301: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1302: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1303: PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);
1304: PetscOptionsGetReal(NULL,NULL,"-sinker_c0x",&opts_c0x,NULL);
1305: PetscOptionsGetReal(NULL,NULL,"-sinker_c0y",&opts_c0y,NULL);
1306: PetscOptionsGetReal(NULL,NULL,"-sinker_s0x",&opts_s0x,NULL);
1307: PetscOptionsGetReal(NULL,NULL,"-sinker_s0y",&opts_s0y,NULL);
1308: PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1309: PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
1310: PetscOptionsGetReal(NULL,NULL,"-sinker_phi",&opts_phi,NULL);
1311: opts_phi *= PETSC_PI / 180;
1313: for (p = 0; p < GAUSS_POINTS; p++) {
1314: coord_x = centroid_x;
1315: coord_y = centroid_y;
1316: if (use_gp_coords) {
1317: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1318: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1319: }
1321: element_props[j][i].eta[p] = opts_eta0;
1322: element_props[j][i].fx[p] = 0.0;
1323: element_props[j][i].fy[p] = -0.2;
1325: radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1326: if (radius2 < opts_r*opts_r
1327: || (PetscAbs(+(coord_x - opts_s0x)*PetscCosReal(opts_phi) + (coord_y - opts_s0y)*PetscSinReal(opts_phi)) < opts_dx/2
1328: && PetscAbs(-(coord_x - opts_s0x)*PetscSinReal(opts_phi) + (coord_y - opts_s0y)*PetscCosReal(opts_phi)) < opts_dy/2)) {
1329: element_props[j][i].eta[p] = opts_eta1;
1330: element_props[j][i].fx[p] = 0.0;
1331: element_props[j][i].fy[p] = -1.0;
1332: }
1333: }
1334: } else if (coefficient_structure == 4) { /* subdomain jump */
1335: PetscReal opts_mag,opts_eta0;
1336: PetscInt opts_nz,px,py;
1337: PetscBool jump;
1339: opts_mag = 1.0;
1340: opts_eta0 = 1.0;
1341: opts_nz = 1;
1343: PetscOptionsGetReal(NULL,NULL,"-jump_eta0",&opts_eta0,NULL);
1344: PetscOptionsGetReal(NULL,NULL,"-jump_magnitude",&opts_mag,NULL);
1345: PetscOptionsGetInt(NULL,NULL,"-jump_nz",&opts_nz,NULL);
1346: DMDAGetInfo(da_Stokes,NULL,NULL,NULL,NULL,&px,&py,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1347: if (px%2) {
1348: jump = (PetscBool)(PetscGlobalRank%2);
1349: } else {
1350: jump = (PetscBool)((PetscGlobalRank/px)%2 ? PetscGlobalRank%2 : !(PetscGlobalRank%2));
1351: }
1352: for (p = 0; p < GAUSS_POINTS; p++) {
1353: coord_x = centroid_x;
1354: coord_y = centroid_y;
1355: if (use_gp_coords) {
1356: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1357: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1358: }
1360: element_props[j][i].eta[p] = jump ? PetscPowReal(10.0,opts_mag) : opts_eta0;
1361: element_props[j][i].fx[p] = 0.0;
1362: element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1363: }
1364: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1365: }
1366: }
1367: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
1369: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1371: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1372: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1373: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
1375: if (output_gnuplot) {
1376: DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1377: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coefficients for Stokes eqn.","properties");
1378: }
1380: if (glvis) {
1381: Vec glv_prop,etaf;
1382: PetscViewer view;
1383: PetscInt dim;
1384: const char *fec = {"FiniteElementCollection: L2_2D_P1"};
1386: DMGetDimension(da_Stokes,&dim);
1387: VecCreateSeq(PETSC_COMM_SELF,GAUSS_POINTS*mx*mx,&etaf);
1388: PetscObjectSetName((PetscObject)etaf,"viscosity");
1389: PetscViewerGLVisOpen(PETSC_COMM_WORLD,PETSC_VIEWER_GLVIS_SOCKET,NULL,PETSC_DECIDE,&view);
1390: PetscViewerGLVisSetFields(view,1,&fec,&dim,glvis_extract_eta,(PetscObject*)&etaf,da_prop,NULL);
1391: DMGetLocalVector(da_prop,&glv_prop);
1392: DMGlobalToLocalBegin(da_prop,properties,INSERT_VALUES,glv_prop);
1393: DMGlobalToLocalEnd(da_prop,properties,INSERT_VALUES,glv_prop);
1394: VecSetDM(etaf,da_Stokes);
1395: VecView(glv_prop,view);
1396: PetscViewerDestroy(&view);
1397: DMRestoreLocalVector(da_prop,&glv_prop);
1398: VecDestroy(&etaf);
1399: }
1401: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1402: DMCreateMatrix(da_Stokes,&A);
1403: DMCreateMatrix(da_Stokes,&B);
1404: DMCreateGlobalVector(da_Stokes,&f);
1405: DMCreateGlobalVector(da_Stokes,&X);
1407: /* assemble A11 */
1408: MatZeroEntries(A);
1409: MatZeroEntries(B);
1410: VecZeroEntries(f);
1412: AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1413: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1414: /* build force vector */
1415: AssembleF_Stokes(f,da_Stokes,da_prop,properties);
1417: DMDABCApplyFreeSlip(da_Stokes,A,f);
1418: DMDABCApplyFreeSlip(da_Stokes,B,NULL);
1420: /* SOLVE */
1421: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1422: KSPSetOptionsPrefix(ksp_S,"stokes_");
1423: KSPSetDM(ksp_S,da_Stokes);
1424: KSPSetDMActive(ksp_S,PETSC_FALSE);
1425: KSPSetOperators(ksp_S,A,B);
1426: KSPSetFromOptions(ksp_S);
1427: {
1428: PC pc;
1429: const PetscInt ufields[] = {0,1},pfields[1] = {2};
1431: KSPGetPC(ksp_S,&pc);
1432: PCFieldSplitSetBlockSize(pc,3);
1433: PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1434: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1435: }
1437: {
1438: PC pc;
1439: PetscBool same = PETSC_FALSE;
1440: KSPGetPC(ksp_S,&pc);
1441: PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&same);
1442: if (same) {
1443: PetscBool usedivmat = PETSC_FALSE;
1444: KSPSetOperators(ksp_S,A,A);
1446: PetscOptionsGetBool(NULL,NULL,"-stokes_pc_bddc_use_divergence",&usedivmat,NULL);
1447: if (usedivmat) {
1448: IS *fields,vel;
1449: PetscInt i,nf;
1451: DMCreateFieldDecomposition(da_Stokes,&nf,NULL,&fields,NULL);
1452: ISConcatenate(PETSC_COMM_WORLD,2,fields,&vel);
1453: MatZeroRowsIS(B,fields[2],1.0,NULL,NULL); /* we put 1.0 on the diagonal to pick the pressure average too */
1454: MatTranspose(B,MAT_INPLACE_MATRIX,&B);
1455: MatZeroRowsIS(B,vel,0.0,NULL,NULL);
1456: ISDestroy(&vel);
1457: PCBDDCSetDivergenceMat(pc,B,PETSC_FALSE,NULL);
1458: for (i=0;i<nf;i++) {
1459: ISDestroy(&fields[i]);
1460: }
1461: PetscFree(fields);
1462: }
1463: }
1464: }
1466: {
1467: PC pc_S;
1468: KSP *sub_ksp,ksp_U;
1469: PetscInt nsplits;
1470: DM da_U;
1471: PetscBool is_pcfs;
1473: KSPSetUp(ksp_S);
1474: KSPGetPC(ksp_S,&pc_S);
1476: is_pcfs = PETSC_FALSE;
1477: PetscObjectTypeCompare((PetscObject)pc_S,PCFIELDSPLIT,&is_pcfs);
1479: if (is_pcfs) {
1480: PCFieldSplitGetSubKSP(pc_S,&nsplits,&sub_ksp);
1481: ksp_U = sub_ksp[0];
1482: PetscFree(sub_ksp);
1484: if (nsplits == 2) {
1485: DMDACreateCompatibleDMDA(da_Stokes,2,&da_U);
1487: KSPSetDM(ksp_U,da_U);
1488: KSPSetDMActive(ksp_U,PETSC_FALSE);
1489: DMDestroy(&da_U);
1490: if (change) {
1491: const char opt[] = "-stokes_fieldsplit_u_pc_factor_mat_solver_type mumps";
1492: PetscOptionsInsertString(NULL,opt);
1493: KSPSetFromOptions(ksp_U);
1494: }
1495: KSPSetUp(ksp_U);
1496: }
1497: }
1498: }
1500: KSPSolve(ksp_S,f,X);
1502: PetscOptionsGetString(NULL,NULL,"-o",filename,sizeof(filename),&set);
1503: if (set) {
1504: char *ext;
1505: PetscViewer viewer;
1506: PetscBool flg;
1507: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1508: PetscStrrchr(filename,'.',&ext);
1509: PetscStrcmp("vts",ext,&flg);
1510: if (flg) {
1511: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1512: } else {
1513: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1514: }
1515: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1516: PetscViewerFileSetName(viewer,filename);
1517: VecView(X,viewer);
1518: PetscViewerDestroy(&viewer);
1519: }
1520: if (output_gnuplot) {
1521: DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
1522: }
1524: if (glvis) {
1525: PetscViewer view;
1527: PetscViewerCreate(PETSC_COMM_WORLD,&view);
1528: PetscViewerSetType(view,PETSCVIEWERGLVIS);
1529: VecView(X,view);
1530: PetscViewerDestroy(&view);
1531: }
1533: KSPGetIterationNumber(ksp_S,&its);
1535: if (coefficient_structure == 0) {
1536: PetscReal opts_eta0,opts_eta1,opts_xc;
1537: PetscInt opts_nz,N;
1538: DM da_Stokes_analytic;
1539: Vec X_analytic;
1540: PetscReal nrm1[3],nrm2[3],nrmI[3];
1542: opts_eta0 = 1.0;
1543: opts_eta1 = 1.0;
1544: opts_xc = 0.5;
1545: opts_nz = 1;
1547: PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1548: PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1549: PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1550: PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);
1552: DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1553: if (output_gnuplot) {
1554: DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
1555: }
1556: DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);
1558: VecAXPY(X_analytic,-1.0,X);
1559: VecGetSize(X_analytic,&N);
1560: N = N/3;
1562: VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1563: VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1564: VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);
1566: VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1567: VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1568: VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);
1570: VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1571: VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1572: VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);
1574: DMDestroy(&da_Stokes_analytic);
1575: VecDestroy(&X_analytic);
1576: }
1578: KSPDestroy(&ksp_S);
1579: VecDestroy(&X);
1580: VecDestroy(&f);
1581: MatDestroy(&A);
1582: MatDestroy(&B);
1584: DMDestroy(&da_Stokes);
1585: DMDestroy(&da_prop);
1587: VecDestroy(&properties);
1588: VecDestroy(&l_properties);
1589: return 0;
1590: }
1592: int main(int argc,char **args)
1593: {
1594: PetscInt mx,my;
1596: PetscInitialize(&argc,&args,(char*)0,help);
1597: mx = my = 10;
1598: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1599: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1600: solve_stokes_2d_coupled(mx,my);
1601: PetscFinalize();
1602: return 0;
1603: }
1605: /* -------------------------- helpers for boundary conditions -------------------------------- */
1606: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
1607: {
1608: DM cda;
1609: Vec coords;
1610: PetscInt si,sj,nx,ny,i,j;
1611: PetscInt M,N;
1612: DMDACoor2d **_coords;
1613: const PetscInt *g_idx;
1614: PetscInt *bc_global_ids;
1615: PetscScalar *bc_vals;
1616: PetscInt nbcs;
1617: PetscInt n_dofs;
1618: ISLocalToGlobalMapping ltogm;
1621: DMGetLocalToGlobalMapping(da,<ogm);
1622: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1624: DMGetCoordinateDM(da,&cda);
1625: DMGetCoordinatesLocal(da,&coords);
1626: DMDAVecGetArray(cda,coords,&_coords);
1627: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1628: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1630: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1631: PetscMalloc1(ny*n_dofs,&bc_vals);
1633: /* init the entries to -1 so VecSetValues will ignore them */
1634: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1636: i = nx-1;
1637: for (j = 0; j < ny; j++) {
1638: PetscInt local_id;
1640: local_id = i+j*nx;
1642: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1644: bc_vals[j] = 0.0;
1645: }
1646: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1647: nbcs = 0;
1648: if ((si+nx) == (M)) nbcs = ny;
1650: if (b) {
1651: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1652: VecAssemblyBegin(b);
1653: VecAssemblyEnd(b);
1654: }
1655: if (A) {
1656: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1657: }
1659: PetscFree(bc_vals);
1660: PetscFree(bc_global_ids);
1662: DMDAVecRestoreArray(cda,coords,&_coords);
1663: return 0;
1664: }
1666: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
1667: {
1668: DM cda;
1669: Vec coords;
1670: PetscInt si,sj,nx,ny,i,j;
1671: PetscInt M,N;
1672: DMDACoor2d **_coords;
1673: const PetscInt *g_idx;
1674: PetscInt *bc_global_ids;
1675: PetscScalar *bc_vals;
1676: PetscInt nbcs;
1677: PetscInt n_dofs;
1678: ISLocalToGlobalMapping ltogm;
1681: DMGetLocalToGlobalMapping(da,<ogm);
1682: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1684: DMGetCoordinateDM(da,&cda);
1685: DMGetCoordinatesLocal(da,&coords);
1686: DMDAVecGetArray(cda,coords,&_coords);
1687: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1688: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1690: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1691: PetscMalloc1(ny*n_dofs,&bc_vals);
1693: /* init the entries to -1 so VecSetValues will ignore them */
1694: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1696: i = 0;
1697: for (j = 0; j < ny; j++) {
1698: PetscInt local_id;
1700: local_id = i+j*nx;
1702: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1704: bc_vals[j] = 0.0;
1705: }
1706: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1707: nbcs = 0;
1708: if (si == 0) nbcs = ny;
1710: if (b) {
1711: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1712: VecAssemblyBegin(b);
1713: VecAssemblyEnd(b);
1714: }
1716: if (A) {
1717: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1718: }
1720: PetscFree(bc_vals);
1721: PetscFree(bc_global_ids);
1723: DMDAVecRestoreArray(cda,coords,&_coords);
1724: return 0;
1725: }
1727: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
1728: {
1729: DM cda;
1730: Vec coords;
1731: PetscInt si,sj,nx,ny,i,j;
1732: PetscInt M,N;
1733: DMDACoor2d **_coords;
1734: const PetscInt *g_idx;
1735: PetscInt *bc_global_ids;
1736: PetscScalar *bc_vals;
1737: PetscInt nbcs;
1738: PetscInt n_dofs;
1739: ISLocalToGlobalMapping ltogm;
1742: DMGetLocalToGlobalMapping(da,<ogm);
1743: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1745: DMGetCoordinateDM(da,&cda);
1746: DMGetCoordinatesLocal(da,&coords);
1747: DMDAVecGetArray(cda,coords,&_coords);
1748: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1749: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1751: PetscMalloc1(nx,&bc_global_ids);
1752: PetscMalloc1(nx,&bc_vals);
1754: /* init the entries to -1 so VecSetValues will ignore them */
1755: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1757: j = ny-1;
1758: for (i = 0; i < nx; i++) {
1759: PetscInt local_id;
1761: local_id = i+j*nx;
1763: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1765: bc_vals[i] = 0.0;
1766: }
1767: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1768: nbcs = 0;
1769: if ((sj+ny) == (N)) nbcs = nx;
1771: if (b) {
1772: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1773: VecAssemblyBegin(b);
1774: VecAssemblyEnd(b);
1775: }
1776: if (A) {
1777: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1778: }
1780: PetscFree(bc_vals);
1781: PetscFree(bc_global_ids);
1783: DMDAVecRestoreArray(cda,coords,&_coords);
1784: return 0;
1785: }
1787: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
1788: {
1789: DM cda;
1790: Vec coords;
1791: PetscInt si,sj,nx,ny,i,j;
1792: PetscInt M,N;
1793: DMDACoor2d **_coords;
1794: const PetscInt *g_idx;
1795: PetscInt *bc_global_ids;
1796: PetscScalar *bc_vals;
1797: PetscInt nbcs;
1798: PetscInt n_dofs;
1799: ISLocalToGlobalMapping ltogm;
1802: DMGetLocalToGlobalMapping(da,<ogm);
1803: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1805: DMGetCoordinateDM(da,&cda);
1806: DMGetCoordinatesLocal(da,&coords);
1807: DMDAVecGetArray(cda,coords,&_coords);
1808: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1809: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1811: PetscMalloc1(nx,&bc_global_ids);
1812: PetscMalloc1(nx,&bc_vals);
1814: /* init the entries to -1 so VecSetValues will ignore them */
1815: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1817: j = 0;
1818: for (i = 0; i < nx; i++) {
1819: PetscInt local_id;
1821: local_id = i+j*nx;
1823: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1825: bc_vals[i] = 0.0;
1826: }
1827: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1828: nbcs = 0;
1829: if (sj == 0) nbcs = nx;
1831: if (b) {
1832: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1833: VecAssemblyBegin(b);
1834: VecAssemblyEnd(b);
1835: }
1836: if (A) {
1837: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1838: }
1840: PetscFree(bc_vals);
1841: PetscFree(bc_global_ids);
1843: DMDAVecRestoreArray(cda,coords,&_coords);
1844: return 0;
1845: }
1847: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
1848: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1849: {
1851: BCApplyZero_NORTH(da_Stokes,1,A,f);
1852: BCApplyZero_EAST(da_Stokes,0,A,f);
1853: BCApplyZero_SOUTH(da_Stokes,1,A,f);
1854: BCApplyZero_WEST(da_Stokes,0,A,f);
1855: return 0;
1856: }
1858: /*TEST
1860: build:
1861: requires: !complex !single
1863: test:
1864: args: -stokes_pc_use_amat -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type lu -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short
1866: testset:
1867: args: -stokes_pc_use_amat -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_fieldsplit_u_ksp_type preonly -stokes_fieldsplit_u_pc_type lu -stokes_fieldsplit_p_ksp_type preonly -stokes_fieldsplit_p_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short
1868: test:
1869: suffix: 2
1870: args:
1871: output_file: output/ex43_1.out
1872: test:
1873: requires: mumps
1874: suffix: 2_mumps
1875: args: -change true -stokes_ksp_view
1876: output_file: output/ex43_2_mumps.out
1877: # mumps INFO,INFOG,RINFO,RINFOG may vary on different archs, so keep just a stable one
1878: filter: egrep -v "(INFOG\([^7]|INFO\(|\[0\])"
1880: test:
1881: suffix: 3
1882: nsize: 4
1883: args: -stokes_pc_use_amat -stokes_ksp_type gcr -stokes_ksp_gcr_restart 60 -stokes_ksp_norm_type unpreconditioned -stokes_ksp_rtol 1.e-2 -c_str 3 -sinker_eta0 1.0 -sinker_eta1 100 -sinker_dx 0.4 -sinker_dy 0.3 -mx 128 -my 128 -stokes_ksp_monitor_short -stokes_pc_type mg -stokes_mg_levels_pc_type fieldsplit -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_levels_pc_fieldsplit_block_size 3 -stokes_mg_levels_pc_fieldsplit_0_fields 0,1 -stokes_mg_levels_pc_fieldsplit_1_fields 2 -stokes_mg_levels_fieldsplit_0_pc_type sor -stokes_mg_levels_fieldsplit_1_pc_type sor -stokes_mg_levels_ksp_type chebyshev -stokes_mg_levels_ksp_max_it 1 -stokes_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -stokes_pc_mg_levels 4 -stokes_ksp_view
1885: test:
1886: suffix: 4
1887: nsize: 4
1888: args: -stokes_ksp_type pipegcr -stokes_ksp_pipegcr_mmax 60 -stokes_ksp_pipegcr_unroll_w 1 -stokes_ksp_norm_type natural -c_str 3 -sinker_eta0 1.0 -sinker_eta1 100 -sinker_dx 0.4 -sinker_dy 0.3 -mx 128 -my 128 -stokes_ksp_monitor_short -stokes_pc_type mg -stokes_mg_levels_pc_type fieldsplit -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_levels_pc_fieldsplit_block_size 3 -stokes_mg_levels_pc_fieldsplit_0_fields 0,1 -stokes_mg_levels_pc_fieldsplit_1_fields 2 -stokes_mg_levels_fieldsplit_0_pc_type sor -stokes_mg_levels_fieldsplit_1_pc_type sor -stokes_mg_levels_ksp_type chebyshev -stokes_mg_levels_ksp_max_it 1 -stokes_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -stokes_pc_mg_levels 4 -stokes_ksp_view
1890: test:
1891: suffix: 5
1892: nsize: 4
1893: args: -stokes_pc_fieldsplit_off_diag_use_amat -stokes_ksp_type pipegcr -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type bjacobi -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type bjacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short -stokes_ksp_view
1895: test:
1896: suffix: 6
1897: nsize: 8
1898: args: -stokes_ksp_view -stokes_pc_type mg -stokes_pc_mg_levels 2 -stokes_mg_coarse_pc_type telescope -stokes_mg_coarse_pc_telescope_reduction_factor 2 -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_coarse_pc_telescope_subcomm_type contiguous
1900: test:
1901: suffix: bjacobi
1902: nsize: 4
1903: args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type aij -stokes_ksp_converged_reason
1905: test:
1906: suffix: bjacobi_baij
1907: nsize: 4
1908: args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type baij -stokes_ksp_converged_reason
1909: output_file: output/ex43_bjacobi.out
1911: test:
1912: suffix: nested_gmg
1913: nsize: 4
1914: args: -stokes_pc_fieldsplit_off_diag_use_amat -mx 16 -my 16 -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_fieldsplit_u_pc_type mg -stokes_fieldsplit_u_pc_mg_levels 5 -stokes_fieldsplit_u_pc_mg_galerkin pmat -stokes_fieldsplit_u_ksp_type cg -stokes_fieldsplit_u_ksp_rtol 1.0e-4 -stokes_fieldsplit_u_mg_levels_pc_type jacobi -solcx_eta0 1.0e4 -stokes_fieldsplit_u_ksp_converged_reason -stokes_ksp_converged_reason -stokes_fieldsplit_p_sub_pc_factor_zeropivot 1.e-8
1916: test:
1917: suffix: fetidp
1918: nsize: 8
1919: args: -dm_mat_type is -stokes_ksp_type fetidp -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_converged_reason -stokes_fetidp_pc_fieldsplit_schur_fact_type diag -stokes_fetidp_fieldsplit_p_pc_type bjacobi -stokes_fetidp_fieldsplit_lag_ksp_type preonly -stokes_fetidp_fieldsplit_p_ksp_type preonly -stokes_ksp_fetidp_pressure_field 2 -stokes_fetidp_pc_fieldsplit_schur_scale -1
1921: test:
1922: suffix: fetidp_unsym
1923: nsize: 8
1924: args: -dm_mat_type is -stokes_ksp_type fetidp -stokes_ksp_monitor_true_residual -stokes_ksp_converged_reason -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
1926: test:
1927: suffix: bddc_stokes_deluxe
1928: nsize: 8
1929: args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_coarse_redundant_pc_type svd -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc
1931: test:
1932: suffix: bddc_stokes_subdomainjump_deluxe
1933: nsize: 9
1934: args: -c_str 4 -jump_magnitude 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_coarse_redundant_pc_type svd -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc
1936: TEST*/