Actual source code: ex43.c
petsc-3.11.4 2019-09-28
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);
64: #define NSD 2 /* number of spatial dimensions */
65: #define NODES_PER_EL 4 /* nodes per element */
66: #define U_DOFS 2 /* degrees of freedom per velocity node */
67: #define P_DOFS 1 /* degrees of freedom per pressure node */
68: #define GAUSS_POINTS 4
70: /* Gauss point based evaluation 8+4+4+4 = 20 */
71: typedef struct {
72: PetscScalar gp_coords[2*GAUSS_POINTS];
73: PetscScalar eta[GAUSS_POINTS];
74: PetscScalar fx[GAUSS_POINTS];
75: PetscScalar fy[GAUSS_POINTS];
76: } GaussPointCoefficients;
78: typedef struct {
79: PetscScalar u_dof;
80: PetscScalar v_dof;
81: PetscScalar p_dof;
82: } StokesDOF;
84: static PetscErrorCode glvis_extract_eta(PetscObject oV,PetscInt nf, PetscObject oVf[], void *ctx)
85: {
86: DM properties_da = (DM)(ctx),stokes_da;
87: Vec V = (Vec)oV, *Vf = (Vec*)oVf;
88: GaussPointCoefficients **props;
89: PetscInt sex,sey,mx,my;
90: PetscInt ei,ej,p,cum;
91: PetscScalar *array;
92: PetscErrorCode ierr;
95: VecGetDM(Vf[0],&stokes_da);
96: DMDAVecGetArrayRead(properties_da,V,&props);
97: DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
98: DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
99: VecGetArray(Vf[0],&array);
100: cum = 0;
101: for (ej = sey; ej < sey+my; ej++) {
102: for (ei = sex; ei < sex+mx; ei++) {
103: for (p = 0; p < GAUSS_POINTS; p++) {
104: array[cum++] = props[ej][ei].eta[p];
105: }
106: }
107: }
108: VecRestoreArray(Vf[0],&array);
109: DMDAVecRestoreArrayRead(properties_da,V,&props);
110: return(0);
111: }
113: /*
114: Element: Local basis function ordering
115: 1-----2
116: | |
117: | |
118: 0-----3
119: */
120: static void ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])
121: {
122: PetscScalar xi = _xi[0];
123: PetscScalar eta = _xi[1];
125: Ni[0] = 0.25*(1.0-xi)*(1.0-eta);
126: Ni[1] = 0.25*(1.0-xi)*(1.0+eta);
127: Ni[2] = 0.25*(1.0+xi)*(1.0+eta);
128: Ni[3] = 0.25*(1.0+xi)*(1.0-eta);
129: }
131: static void ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
132: {
133: PetscScalar xi = _xi[0];
134: PetscScalar eta = _xi[1];
136: GNi[0][0] = -0.25*(1.0-eta);
137: GNi[0][1] = -0.25*(1.0+eta);
138: GNi[0][2] = 0.25*(1.0+eta);
139: GNi[0][3] = 0.25*(1.0-eta);
141: GNi[1][0] = -0.25*(1.0-xi);
142: GNi[1][1] = 0.25*(1.0-xi);
143: GNi[1][2] = 0.25*(1.0+xi);
144: GNi[1][3] = -0.25*(1.0+xi);
145: }
147: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
148: {
149: PetscScalar J00,J01,J10,J11,J;
150: PetscScalar iJ00,iJ01,iJ10,iJ11;
151: PetscInt i;
153: J00 = J01 = J10 = J11 = 0.0;
154: for (i = 0; i < NODES_PER_EL; i++) {
155: PetscScalar cx = coords[2*i];
156: PetscScalar cy = coords[2*i+1];
158: J00 = J00+GNi[0][i]*cx; /* J_xx = dx/dxi */
159: J01 = J01+GNi[0][i]*cy; /* J_xy = dy/dxi */
160: J10 = J10+GNi[1][i]*cx; /* J_yx = dx/deta */
161: J11 = J11+GNi[1][i]*cy; /* J_yy = dy/deta */
162: }
163: J = (J00*J11)-(J01*J10);
165: iJ00 = J11/J;
166: iJ01 = -J01/J;
167: iJ10 = -J10/J;
168: iJ11 = J00/J;
170: for (i = 0; i < NODES_PER_EL; i++) {
171: GNx[0][i] = GNi[0][i]*iJ00+GNi[1][i]*iJ01;
172: GNx[1][i] = GNi[0][i]*iJ10+GNi[1][i]*iJ11;
173: }
175: if (det_J) *det_J = J;
176: }
178: static void ConstructGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
179: {
180: *ngp = 4;
181: gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919;
182: gp_xi[1][0] = -0.57735026919; gp_xi[1][1] = 0.57735026919;
183: gp_xi[2][0] = 0.57735026919; gp_xi[2][1] = 0.57735026919;
184: gp_xi[3][0] = 0.57735026919; gp_xi[3][1] = -0.57735026919;
185: gp_weight[0] = 1.0;
186: gp_weight[1] = 1.0;
187: gp_weight[2] = 1.0;
188: gp_weight[3] = 1.0;
189: }
191: /*
192: i,j are the element indices
193: The unknown is a vector quantity.
194: The s[].c is used to indicate the degree of freedom.
195: */
196: static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)
197: {
199: /* velocity */
200: /* node 0 */
201: s_u[0].i = i; s_u[0].j = j; s_u[0].c = 0; /* Vx0 */
202: s_u[1].i = i; s_u[1].j = j; s_u[1].c = 1; /* Vy0 */
204: /* node 1 */
205: s_u[2].i = i; s_u[2].j = j+1; s_u[2].c = 0; /* Vx1 */
206: s_u[3].i = i; s_u[3].j = j+1; s_u[3].c = 1; /* Vy1 */
208: /* node 2 */
209: s_u[4].i = i+1; s_u[4].j = j+1; s_u[4].c = 0; /* Vx2 */
210: s_u[5].i = i+1; s_u[5].j = j+1; s_u[5].c = 1; /* Vy2 */
212: /* node 3 */
213: s_u[6].i = i+1; s_u[6].j = j; s_u[6].c = 0; /* Vx3 */
214: s_u[7].i = i+1; s_u[7].j = j; s_u[7].c = 1; /* Vy3 */
216: /* pressure */
217: s_p[0].i = i; s_p[0].j = j; s_p[0].c = 2; /* P0 */
218: s_p[1].i = i; s_p[1].j = j+1; s_p[1].c = 2; /* P1 */
219: s_p[2].i = i+1; s_p[2].j = j+1; s_p[2].c = 2; /* P2 */
220: s_p[3].i = i+1; s_p[3].j = j; s_p[3].c = 2; /* P3 */
221: return(0);
222: }
224: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da,PetscInt **_lx,PetscInt **_ly)
225: {
227: PetscMPIInt rank;
228: PetscInt proc_I,proc_J;
229: PetscInt cpu_x,cpu_y;
230: PetscInt local_mx,local_my;
231: Vec vlx,vly;
232: PetscInt *LX,*LY,i;
233: PetscScalar *_a;
234: Vec V_SEQ;
235: VecScatter ctx;
238: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
240: DMDAGetInfo(da,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
242: proc_J = rank/cpu_x;
243: proc_I = rank-cpu_x*proc_J;
245: PetscMalloc1(cpu_x,&LX);
246: PetscMalloc1(cpu_y,&LY);
248: DMDAGetElementsSizes(da,&local_mx,&local_my,NULL);
249: VecCreate(PETSC_COMM_WORLD,&vlx);
250: VecSetSizes(vlx,PETSC_DECIDE,cpu_x);
251: VecSetFromOptions(vlx);
253: VecCreate(PETSC_COMM_WORLD,&vly);
254: VecSetSizes(vly,PETSC_DECIDE,cpu_y);
255: VecSetFromOptions(vly);
257: VecSetValue(vlx,proc_I,(PetscScalar)(local_mx+1.0e-9),INSERT_VALUES);
258: VecSetValue(vly,proc_J,(PetscScalar)(local_my+1.0e-9),INSERT_VALUES);
259: VecAssemblyBegin(vlx);VecAssemblyEnd(vlx);
260: VecAssemblyBegin(vly);VecAssemblyEnd(vly);
262: VecScatterCreateToAll(vlx,&ctx,&V_SEQ);
263: VecScatterBegin(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
264: VecScatterEnd(ctx,vlx,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
265: VecGetArray(V_SEQ,&_a);
266: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
267: VecRestoreArray(V_SEQ,&_a);
268: VecScatterDestroy(&ctx);
269: VecDestroy(&V_SEQ);
271: VecScatterCreateToAll(vly,&ctx,&V_SEQ);
272: VecScatterBegin(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
273: VecScatterEnd(ctx,vly,V_SEQ,INSERT_VALUES,SCATTER_FORWARD);
274: VecGetArray(V_SEQ,&_a);
275: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
276: VecRestoreArray(V_SEQ,&_a);
277: VecScatterDestroy(&ctx);
278: VecDestroy(&V_SEQ);
280: *_lx = LX;
281: *_ly = LY;
283: VecDestroy(&vlx);
284: VecDestroy(&vly);
285: return(0);
286: }
288: static PetscErrorCode DMDACoordViewGnuplot2d(DM da,const char prefix[])
289: {
290: DM cda;
291: Vec coords;
292: DMDACoor2d **_coords;
293: PetscInt si,sj,nx,ny,i,j;
294: FILE *fp;
295: char fname[PETSC_MAX_PATH_LEN];
296: PetscMPIInt rank;
300: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
301: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
302: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
303: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
305: PetscFPrintf(PETSC_COMM_SELF,fp,"### Element geometry for processor %1.4d ### \n",rank);
307: DMGetCoordinateDM(da,&cda);
308: DMGetCoordinatesLocal(da,&coords);
309: DMDAVecGetArrayRead(cda,coords,&_coords);
310: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
311: for (j = sj; j < sj+ny-1; j++) {
312: for (i = si; i < si+nx-1; i++) {
313: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
314: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j+1][i].x),(double)PetscRealPart(_coords[j+1][i].y));
315: 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));
316: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n",(double)PetscRealPart(_coords[j][i+1].x),(double)PetscRealPart(_coords[j][i+1].y));
317: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e \n\n",(double)PetscRealPart(_coords[j][i].x),(double)PetscRealPart(_coords[j][i].y));
318: }
319: }
320: DMDAVecRestoreArrayRead(cda,coords,&_coords);
322: PetscFClose(PETSC_COMM_SELF,fp);
323: return(0);
324: }
326: static PetscErrorCode DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
327: {
328: DM cda;
329: Vec coords,local_fields;
330: DMDACoor2d **_coords;
331: FILE *fp;
332: char fname[PETSC_MAX_PATH_LEN];
333: PetscMPIInt rank;
334: PetscInt si,sj,nx,ny,i,j;
335: PetscInt n_dofs,d;
336: PetscScalar *_fields;
340: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
341: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
342: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
343: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
345: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
346: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
347: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
348: for (d = 0; d < n_dofs; d++) {
349: const char *field_name;
350: DMDAGetFieldName(da,d,&field_name);
351: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
352: }
353: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
355: DMGetCoordinateDM(da,&cda);
356: DMGetCoordinatesLocal(da,&coords);
357: DMDAVecGetArray(cda,coords,&_coords);
358: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
360: DMCreateLocalVector(da,&local_fields);
361: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
362: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
363: VecGetArray(local_fields,&_fields);
365: for (j = sj; j < sj+ny; j++) {
366: for (i = si; i < si+nx; i++) {
367: PetscScalar coord_x,coord_y;
368: PetscScalar field_d;
370: coord_x = _coords[j][i].x;
371: coord_y = _coords[j][i].y;
373: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
374: for (d = 0; d < n_dofs; d++) {
375: field_d = _fields[n_dofs*((i-si)+(j-sj)*(nx))+d];
376: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e ",(double)PetscRealPart(field_d));
377: }
378: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
379: }
380: }
381: VecRestoreArray(local_fields,&_fields);
382: VecDestroy(&local_fields);
384: DMDAVecRestoreArray(cda,coords,&_coords);
386: PetscFClose(PETSC_COMM_SELF,fp);
387: return(0);
388: }
390: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])
391: {
392: DM cda;
393: Vec local_fields;
394: FILE *fp;
395: char fname[PETSC_MAX_PATH_LEN];
396: PetscMPIInt rank;
397: PetscInt si,sj,nx,ny,i,j,p;
398: PetscInt n_dofs,d;
399: GaussPointCoefficients **_coefficients;
400: PetscErrorCode ierr;
403: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
404: PetscSNPrintf(fname,sizeof(fname),"%s-p%1.4d.dat",prefix,rank);
405: PetscFOpen(PETSC_COMM_SELF,fname,"w",&fp);
406: if (!fp) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file");
408: PetscFPrintf(PETSC_COMM_SELF,fp,"### %s (processor %1.4d) ### \n",comment,rank);
409: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_dofs,0,0,0,0,0);
410: PetscFPrintf(PETSC_COMM_SELF,fp,"### x y ");
411: for (d = 0; d < n_dofs; d++) {
412: const char *field_name;
413: DMDAGetFieldName(da,d,&field_name);
414: PetscFPrintf(PETSC_COMM_SELF,fp,"%s ",field_name);
415: }
416: PetscFPrintf(PETSC_COMM_SELF,fp,"###\n");
418: DMGetCoordinateDM(da,&cda);
419: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
421: DMCreateLocalVector(da,&local_fields);
422: DMGlobalToLocalBegin(da,fields,INSERT_VALUES,local_fields);
423: DMGlobalToLocalEnd(da,fields,INSERT_VALUES,local_fields);
424: DMDAVecGetArray(da,local_fields,&_coefficients);
426: for (j = sj; j < sj+ny; j++) {
427: for (i = si; i < si+nx; i++) {
428: PetscScalar coord_x,coord_y;
430: for (p = 0; p < GAUSS_POINTS; p++) {
431: coord_x = _coefficients[j][i].gp_coords[2*p];
432: coord_y = _coefficients[j][i].gp_coords[2*p+1];
434: PetscFPrintf(PETSC_COMM_SELF,fp,"%1.6e %1.6e ",(double)PetscRealPart(coord_x),(double)PetscRealPart(coord_y));
436: 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]));
437: PetscFPrintf(PETSC_COMM_SELF,fp,"\n");
438: }
439: }
440: }
441: DMDAVecRestoreArray(da,local_fields,&_coefficients);
442: VecDestroy(&local_fields);
444: PetscFClose(PETSC_COMM_SELF,fp);
445: return(0);
446: }
448: 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)
449: {
450: PetscInt ij;
451: PetscInt r,c,nc;
453: nc = u_NPE*u_dof;
454: r = w_dof*wi+wd;
455: c = u_dof*ui+ud;
456: ij = r*nc+c;
457: return ij;
458: }
460: /*
461: D = [ 2.eta 0 0 ]
462: [ 0 2.eta 0 ]
463: [ 0 0 eta ]
465: B = [ d_dx 0 ]
466: [ 0 d_dy ]
467: [ d_dy d_dx ]
468: */
469: static void FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
470: {
471: PetscInt ngp;
472: PetscScalar gp_xi[GAUSS_POINTS][2];
473: PetscScalar gp_weight[GAUSS_POINTS];
474: PetscInt p,i,j,k;
475: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
476: PetscScalar J_p,tildeD[3];
477: PetscScalar B[3][U_DOFS*NODES_PER_EL];
479: /* define quadrature rule */
480: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
482: /* evaluate integral */
483: for (p = 0; p < ngp; p++) {
484: ConstructQ12D_GNi(gp_xi[p],GNi_p);
485: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
487: for (i = 0; i < NODES_PER_EL; i++) {
488: PetscScalar d_dx_i = GNx_p[0][i];
489: PetscScalar d_dy_i = GNx_p[1][i];
491: B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
492: B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
493: B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
494: }
496: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
497: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
498: tildeD[2] = gp_weight[p]*J_p*eta[p];
500: /* form Bt tildeD B */
501: /*
502: Ke_ij = Bt_ik . D_kl . B_lj
503: = B_ki . D_kl . B_lj
504: = B_ki . D_kk . B_kj
505: */
506: for (i = 0; i < 8; i++) {
507: for (j = 0; j < 8; j++) {
508: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
509: Ke[i+8*j] = Ke[i+8*j]+B[k][i]*tildeD[k]*B[k][j];
510: }
511: }
512: }
513: }
514: }
516: static void FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])
517: {
518: PetscInt ngp;
519: PetscScalar gp_xi[GAUSS_POINTS][2];
520: PetscScalar gp_weight[GAUSS_POINTS];
521: PetscInt p,i,j,di;
522: PetscScalar Ni_p[NODES_PER_EL];
523: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
524: PetscScalar J_p,fac;
526: /* define quadrature rule */
527: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
529: /* evaluate integral */
530: for (p = 0; p < ngp; p++) {
531: ConstructQ12D_Ni(gp_xi[p],Ni_p);
532: ConstructQ12D_GNi(gp_xi[p],GNi_p);
533: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
534: fac = gp_weight[p]*J_p;
536: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
537: for (di = 0; di < NSD; di++) { /* u dofs */
538: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
539: PetscInt IJ;
540: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);
542: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
543: }
544: }
545: }
546: }
547: }
549: static void FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])
550: {
551: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
552: PetscInt i,j;
553: PetscInt nr_g,nc_g;
555: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
556: FormGradientOperatorQ1(Ge,coords);
558: nr_g = U_DOFS*NODES_PER_EL;
559: nc_g = P_DOFS*NODES_PER_EL;
561: for (i = 0; i < nr_g; i++) {
562: for (j = 0; j < nc_g; j++) {
563: De[nr_g*j+i] = Ge[nc_g*i+j];
564: }
565: }
566: }
568: static void FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
569: {
570: PetscInt ngp;
571: PetscScalar gp_xi[GAUSS_POINTS][2];
572: PetscScalar gp_weight[GAUSS_POINTS];
573: PetscInt p,i,j;
574: PetscScalar Ni_p[NODES_PER_EL];
575: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
576: PetscScalar J_p,fac,eta_avg;
578: /* define quadrature rule */
579: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
581: /* evaluate integral */
582: for (p = 0; p < ngp; p++) {
583: ConstructQ12D_Ni(gp_xi[p],Ni_p);
584: ConstructQ12D_GNi(gp_xi[p],GNi_p);
585: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
586: fac = gp_weight[p]*J_p;
588: for (i = 0; i < NODES_PER_EL; i++) {
589: for (j = 0; j < NODES_PER_EL; j++) {
590: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*(Ni_p[i]*Ni_p[j]-0.0625);
591: }
592: }
593: }
595: /* scale */
596: eta_avg = 0.0;
597: for (p = 0; p < ngp; p++) eta_avg += eta[p];
598: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
599: fac = 1.0/eta_avg;
600: for (i = 0; i < NODES_PER_EL; i++) {
601: for (j = 0; j < NODES_PER_EL; j++) {
602: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
603: }
604: }
605: }
607: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
608: {
609: PetscInt ngp;
610: PetscScalar gp_xi[GAUSS_POINTS][2];
611: PetscScalar gp_weight[GAUSS_POINTS];
612: PetscInt p,i,j;
613: PetscScalar Ni_p[NODES_PER_EL];
614: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
615: PetscScalar J_p,fac,eta_avg;
617: /* define quadrature rule */
618: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
620: /* evaluate integral */
621: for (p = 0; p < ngp; p++) {
622: ConstructQ12D_Ni(gp_xi[p],Ni_p);
623: ConstructQ12D_GNi(gp_xi[p],GNi_p);
624: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
625: fac = gp_weight[p]*J_p;
627: for (i = 0; i < NODES_PER_EL; i++) {
628: for (j = 0; j < NODES_PER_EL; j++) {
629: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
630: }
631: }
632: }
634: /* scale */
635: eta_avg = 0.0;
636: for (p = 0; p < ngp; p++) eta_avg += eta[p];
637: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
638: fac = 1.0/eta_avg;
639: for (i = 0; i < NODES_PER_EL; i++) {
640: for (j = 0; j < NODES_PER_EL; j++) {
641: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
642: }
643: }
644: }
646: static void FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
647: {
648: PetscInt ngp;
649: PetscScalar gp_xi[GAUSS_POINTS][2];
650: PetscScalar gp_weight[GAUSS_POINTS];
651: PetscInt p,i;
652: PetscScalar Ni_p[NODES_PER_EL];
653: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
654: PetscScalar J_p,fac;
656: /* define quadrature rule */
657: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
659: /* evaluate integral */
660: for (p = 0; p < ngp; p++) {
661: ConstructQ12D_Ni(gp_xi[p],Ni_p);
662: ConstructQ12D_GNi(gp_xi[p],GNi_p);
663: ConstructQ12D_GNx(GNi_p,GNx_p,coords,&J_p);
664: fac = gp_weight[p]*J_p;
666: for (i = 0; i < NODES_PER_EL; i++) {
667: Fe[NSD*i] += fac*Ni_p[i]*fx[p];
668: Fe[NSD*i+1] += fac*Ni_p[i]*fy[p];
669: }
670: }
671: }
673: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])
674: {
676: /* get coords for the element */
677: el_coords[NSD*0] = _coords[ej][ei].x; el_coords[NSD*0+1] = _coords[ej][ei].y;
678: el_coords[NSD*1] = _coords[ej+1][ei].x; el_coords[NSD*1+1] = _coords[ej+1][ei].y;
679: el_coords[NSD*2] = _coords[ej+1][ei+1].x; el_coords[NSD*2+1] = _coords[ej+1][ei+1].y;
680: el_coords[NSD*3] = _coords[ej][ei+1].x; el_coords[NSD*3+1] = _coords[ej][ei+1].y;
681: return(0);
682: }
684: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
685: {
686: DM cda;
687: Vec coords;
688: DMDACoor2d **_coords;
689: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
690: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
691: PetscInt sex,sey,mx,my;
692: PetscInt ei,ej;
693: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
694: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
695: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
696: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
697: PetscScalar el_coords[NODES_PER_EL*NSD];
698: Vec local_properties;
699: GaussPointCoefficients **props;
700: PetscScalar *prop_eta;
701: PetscErrorCode ierr;
704: /* setup for coords */
705: DMGetCoordinateDM(stokes_da,&cda);
706: DMGetCoordinatesLocal(stokes_da,&coords);
707: DMDAVecGetArray(cda,coords,&_coords);
709: /* setup for coefficients */
710: DMCreateLocalVector(properties_da,&local_properties);
711: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
712: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
713: DMDAVecGetArray(properties_da,local_properties,&props);
715: DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
716: DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
717: for (ej = sey; ej < sey+my; ej++) {
718: for (ei = sex; ei < sex+mx; ei++) {
719: /* get coords for the element */
720: GetElementCoords(_coords,ei,ej,el_coords);
722: /* get coefficients for the element */
723: prop_eta = props[ej][ei].eta;
725: /* initialise element stiffness matrix */
726: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
727: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
728: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
729: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
731: /* form element stiffness matrix */
732: FormStressOperatorQ1(Ae,el_coords,prop_eta);
733: FormGradientOperatorQ1(Ge,el_coords);
734: FormDivergenceOperatorQ1(De,el_coords);
735: FormStabilisationOperatorQ1(Ce,el_coords,prop_eta);
737: /* insert element matrix into global matrix */
738: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
739: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
740: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
741: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
742: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
743: }
744: }
745: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
746: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
748: DMDAVecRestoreArray(cda,coords,&_coords);
750: DMDAVecRestoreArray(properties_da,local_properties,&props);
751: VecDestroy(&local_properties);
752: return(0);
753: }
755: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)
756: {
757: DM cda;
758: Vec coords;
759: DMDACoor2d **_coords;
760: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
761: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
762: PetscInt sex,sey,mx,my;
763: PetscInt ei,ej;
764: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
765: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
766: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
767: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
768: PetscScalar el_coords[NODES_PER_EL*NSD];
769: Vec local_properties;
770: GaussPointCoefficients **props;
771: PetscScalar *prop_eta;
772: PetscErrorCode ierr;
775: /* setup for coords */
776: DMGetCoordinateDM(stokes_da,&cda);
777: DMGetCoordinatesLocal(stokes_da,&coords);
778: DMDAVecGetArray(cda,coords,&_coords);
780: /* setup for coefficients */
781: DMCreateLocalVector(properties_da,&local_properties);
782: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
783: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
784: DMDAVecGetArray(properties_da,local_properties,&props);
786: DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
787: DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
788: for (ej = sey; ej < sey+my; ej++) {
789: for (ei = sex; ei < sex+mx; ei++) {
790: /* get coords for the element */
791: GetElementCoords(_coords,ei,ej,el_coords);
793: /* get coefficients for the element */
794: prop_eta = props[ej][ei].eta;
796: /* initialise element stiffness matrix */
797: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
798: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
799: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
800: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
802: /* form element stiffness matrix */
803: FormStressOperatorQ1(Ae,el_coords,prop_eta);
804: FormGradientOperatorQ1(Ge,el_coords);
805: FormScaledMassMatrixOperatorQ1(Ce,el_coords,prop_eta);
807: /* insert element matrix into global matrix */
808: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
809: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
810: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
811: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
812: }
813: }
814: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
815: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
817: DMDAVecRestoreArray(cda,coords,&_coords);
819: DMDAVecRestoreArray(properties_da,local_properties,&props);
820: VecDestroy(&local_properties);
821: return(0);
822: }
824: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
825: {
826: PetscInt n;
829: for (n = 0; n < 4; n++) {
830: 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];
831: 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];
832: 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];
833: }
834: return(0);
835: }
837: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)
838: {
839: DM cda;
840: Vec coords;
841: DMDACoor2d **_coords;
842: MatStencil u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
843: MatStencil p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
844: PetscInt sex,sey,mx,my;
845: PetscInt ei,ej;
846: PetscScalar Fe[NODES_PER_EL*U_DOFS];
847: PetscScalar He[NODES_PER_EL*P_DOFS];
848: PetscScalar el_coords[NODES_PER_EL*NSD];
849: Vec local_properties;
850: GaussPointCoefficients **props;
851: PetscScalar *prop_fx,*prop_fy;
852: Vec local_F;
853: StokesDOF **ff;
854: PetscErrorCode ierr;
857: /* setup for coords */
858: DMGetCoordinateDM(stokes_da,&cda);
859: DMGetCoordinatesLocal(stokes_da,&coords);
860: DMDAVecGetArray(cda,coords,&_coords);
862: /* setup for coefficients */
863: DMGetLocalVector(properties_da,&local_properties);
864: DMGlobalToLocalBegin(properties_da,properties,INSERT_VALUES,local_properties);
865: DMGlobalToLocalEnd(properties_da,properties,INSERT_VALUES,local_properties);
866: DMDAVecGetArray(properties_da,local_properties,&props);
868: /* get acces to the vector */
869: DMGetLocalVector(stokes_da,&local_F);
870: VecZeroEntries(local_F);
871: DMDAVecGetArray(stokes_da,local_F,&ff);
873: DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
874: DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
875: for (ej = sey; ej < sey+my; ej++) {
876: for (ei = sex; ei < sex+mx; ei++) {
877: /* get coords for the element */
878: GetElementCoords(_coords,ei,ej,el_coords);
880: /* get coefficients for the element */
881: prop_fx = props[ej][ei].fx;
882: prop_fy = props[ej][ei].fy;
884: /* initialise element stiffness matrix */
885: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
886: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
888: /* form element stiffness matrix */
889: FormMomentumRhsQ1(Fe,el_coords,prop_fx,prop_fy);
891: /* insert element matrix into global matrix */
892: DMDAGetElementEqnums_up(u_eqn,p_eqn,ei,ej);
894: DMDASetValuesLocalStencil_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
895: }
896: }
898: DMDAVecRestoreArray(stokes_da,local_F,&ff);
899: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
900: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
901: DMRestoreLocalVector(stokes_da,&local_F);
903: DMDAVecRestoreArray(cda,coords,&_coords);
905: DMDAVecRestoreArray(properties_da,local_properties,&props);
906: DMRestoreLocalVector(properties_da,&local_properties);
907: return(0);
908: }
910: static PetscErrorCode DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM *_da,Vec *_X)
911: {
912: DM da,cda;
913: Vec X;
914: StokesDOF **_stokes;
915: Vec coords;
916: DMDACoor2d **_coords;
917: PetscInt si,sj,ei,ej,i,j;
921: 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);
922: DMSetFromOptions(da);
923: DMSetUp(da);
924: DMDASetFieldName(da,0,"anlytic_Vx");
925: DMDASetFieldName(da,1,"anlytic_Vy");
926: DMDASetFieldName(da,2,"analytic_P");
928: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.,0.);
930: DMGetCoordinatesLocal(da,&coords);
931: DMGetCoordinateDM(da,&cda);
932: DMDAVecGetArray(cda,coords,&_coords);
934: DMCreateGlobalVector(da,&X);
935: DMDAVecGetArray(da,X,&_stokes);
937: DMDAGetCorners(da,&si,&sj,0,&ei,&ej,0);
938: for (j = sj; j < sj+ej; j++) {
939: for (i = si; i < si+ei; i++) {
940: PetscReal pos[2],pressure,vel[2],total_stress[3],strain_rate[3];
942: pos[0] = PetscRealPart(_coords[j][i].x);
943: pos[1] = PetscRealPart(_coords[j][i].y);
945: evaluate_solCx(pos,eta0,eta1,xc,nz,vel,&pressure,total_stress,strain_rate);
947: _stokes[j][i].u_dof = vel[0];
948: _stokes[j][i].v_dof = vel[1];
949: _stokes[j][i].p_dof = pressure;
950: }
951: }
952: DMDAVecRestoreArray(da,X,&_stokes);
953: DMDAVecRestoreArray(cda,coords,&_coords);
955: *_da = da;
956: *_X = X;
957: return(0);
958: }
960: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])
961: {
963: /* get the nodal fields */
964: 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;
965: 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;
966: 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;
967: 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;
968: return(0);
969: }
971: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)
972: {
973: DM cda;
974: Vec coords,X_analytic_local,X_local;
975: DMDACoor2d **_coords;
976: PetscInt sex,sey,mx,my;
977: PetscInt ei,ej;
978: PetscScalar el_coords[NODES_PER_EL*NSD];
979: StokesDOF **stokes_analytic,**stokes;
980: StokesDOF stokes_analytic_e[4],stokes_e[4];
982: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
983: PetscScalar Ni_p[NODES_PER_EL];
984: PetscInt ngp;
985: PetscScalar gp_xi[GAUSS_POINTS][2];
986: PetscScalar gp_weight[GAUSS_POINTS];
987: PetscInt p,i;
988: PetscScalar J_p,fac;
989: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
990: PetscInt M;
991: PetscReal xymin[2],xymax[2];
995: /* define quadrature rule */
996: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
998: /* setup for coords */
999: DMGetCoordinateDM(stokes_da,&cda);
1000: DMGetCoordinatesLocal(stokes_da,&coords);
1001: DMDAVecGetArray(cda,coords,&_coords);
1003: /* setup for analytic */
1004: DMCreateLocalVector(stokes_da,&X_analytic_local);
1005: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1006: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1007: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1009: /* setup for solution */
1010: DMCreateLocalVector(stokes_da,&X_local);
1011: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1012: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1013: DMDAVecGetArray(stokes_da,X_local,&stokes);
1015: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1016: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1018: h = (xymax[0]-xymin[0])/((PetscReal)M);
1020: tp_L2 = tu_L2 = tu_H1 = 0.0;
1022: DMDAGetElementsCorners(stokes_da,&sex,&sey,NULL);
1023: DMDAGetElementsSizes(stokes_da,&mx,&my,NULL);
1024: for (ej = sey; ej < sey+my; ej++) {
1025: for (ei = sex; ei < sex+mx; ei++) {
1026: /* get coords for the element */
1027: GetElementCoords(_coords,ei,ej,el_coords);
1028: StokesDAGetNodalFields(stokes,ei,ej,stokes_e);
1029: StokesDAGetNodalFields(stokes_analytic,ei,ej,stokes_analytic_e);
1031: /* evaluate integral */
1032: p_e_L2 = 0.0;
1033: u_e_L2 = 0.0;
1034: u_e_H1 = 0.0;
1035: for (p = 0; p < ngp; p++) {
1036: ConstructQ12D_Ni(gp_xi[p],Ni_p);
1037: ConstructQ12D_GNi(gp_xi[p],GNi_p);
1038: ConstructQ12D_GNx(GNi_p,GNx_p,el_coords,&J_p);
1039: fac = gp_weight[p]*J_p;
1041: for (i = 0; i < NODES_PER_EL; i++) {
1042: PetscScalar u_error,v_error;
1044: 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);
1046: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1047: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1048: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error);
1050: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1051: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error /* du/dy */
1052: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1053: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error); /* dv/dy */
1054: }
1055: }
1057: tp_L2 += p_e_L2;
1058: tu_L2 += u_e_L2;
1059: tu_H1 += u_e_H1;
1060: }
1061: }
1062: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1063: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1064: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1065: p_L2 = PetscSqrtScalar(p_L2);
1066: u_L2 = PetscSqrtScalar(u_L2);
1067: u_H1 = PetscSqrtScalar(u_H1);
1069: 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));
1071: DMDAVecRestoreArray(cda,coords,&_coords);
1073: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1074: VecDestroy(&X_analytic_local);
1075: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1076: VecDestroy(&X_local);
1077: return(0);
1078: }
1080: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx,PetscInt my)
1081: {
1082: DM da_Stokes,da_prop;
1083: PetscInt u_dof,p_dof,dof,stencil_width;
1084: Mat A,B;
1085: DM prop_cda,vel_cda;
1086: Vec prop_coords,vel_coords;
1087: PetscInt si,sj,nx,ny,i,j,p;
1088: Vec f,X;
1089: PetscInt prop_dof,prop_stencil_width;
1090: Vec properties,l_properties;
1091: PetscReal dx,dy;
1092: PetscInt M,N;
1093: DMDACoor2d **_prop_coords,**_vel_coords;
1094: GaussPointCoefficients **element_props;
1095: PetscInt its;
1096: KSP ksp_S;
1097: PetscInt coefficient_structure = 0;
1098: PetscInt cpu_x,cpu_y,*lx = NULL,*ly = NULL;
1099: PetscBool use_gp_coords = PETSC_FALSE,set,output_gnuplot = PETSC_FALSE,glvis = PETSC_FALSE;
1100: char filename[PETSC_MAX_PATH_LEN];
1101: PetscErrorCode ierr;
1105: PetscOptionsGetBool(NULL,NULL,"-gnuplot",&output_gnuplot,NULL);
1106: PetscOptionsGetBool(NULL,NULL,"-glvis",&glvis,NULL);
1108: /* Generate the da for velocity and pressure */
1109: /*
1110: We use Q1 elements for the temperature.
1111: FEM has a 9-point stencil (BOX) or connectivity pattern
1112: Num nodes in each direction is mx+1, my+1
1113: */
1114: u_dof = U_DOFS; /* Vx, Vy - velocities */
1115: p_dof = P_DOFS; /* p - pressure */
1116: dof = u_dof+p_dof;
1117: stencil_width = 1;
1118: 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);
1120: DMSetMatType(da_Stokes,MATAIJ);
1121: DMSetFromOptions(da_Stokes);
1122: DMSetUp(da_Stokes);
1123: DMDASetFieldName(da_Stokes,0,"Vx");
1124: DMDASetFieldName(da_Stokes,1,"Vy");
1125: DMDASetFieldName(da_Stokes,2,"P");
1127: /* unit box [0,1] x [0,1] */
1128: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.,0.);
1130: /* Generate element properties, we will assume all material properties are constant over the element */
1131: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1132: DMDAGetInfo(da_Stokes,0,0,0,0,&cpu_x,&cpu_y,0,0,0,0,0,0,0);
1133: DMDAGetElementOwnershipRanges2d(da_Stokes,&lx,&ly);
1135: prop_dof = (int)(sizeof(GaussPointCoefficients)/sizeof(PetscScalar)); /* gauss point setup */
1136: prop_stencil_width = 0;
1137: 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);
1138: DMSetFromOptions(da_prop);
1139: DMSetUp(da_prop);
1140: PetscFree(lx);
1141: PetscFree(ly);
1143: /* define centroid positions */
1144: DMDAGetInfo(da_prop,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
1145: dx = 1.0/((PetscReal)(M));
1146: dy = 1.0/((PetscReal)(N));
1148: 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);
1150: /* define coefficients */
1151: PetscOptionsGetInt(NULL,NULL,"-c_str",&coefficient_structure,NULL);
1153: DMCreateGlobalVector(da_prop,&properties);
1154: DMCreateLocalVector(da_prop,&l_properties);
1155: DMDAVecGetArray(da_prop,l_properties,&element_props);
1157: DMGetCoordinateDM(da_prop,&prop_cda);
1158: DMGetCoordinatesLocal(da_prop,&prop_coords);
1159: DMDAVecGetArray(prop_cda,prop_coords,&_prop_coords);
1161: DMDAGetGhostCorners(prop_cda,&si,&sj,0,&nx,&ny,0);
1163: DMGetCoordinateDM(da_Stokes,&vel_cda);
1164: DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1165: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1167: /* interpolate the coordinates */
1168: for (j = sj; j < sj+ny; j++) {
1169: for (i = si; i < si+nx; i++) {
1170: PetscInt ngp;
1171: PetscScalar gp_xi[GAUSS_POINTS][2],gp_weight[GAUSS_POINTS];
1172: PetscScalar el_coords[8];
1174: GetElementCoords(_vel_coords,i,j,el_coords);
1175: ConstructGaussQuadrature(&ngp,gp_xi,gp_weight);
1177: for (p = 0; p < GAUSS_POINTS; p++) {
1178: PetscScalar gp_x,gp_y;
1179: PetscInt n;
1180: PetscScalar xi_p[2],Ni_p[4];
1182: xi_p[0] = gp_xi[p][0];
1183: xi_p[1] = gp_xi[p][1];
1184: ConstructQ12D_Ni(xi_p,Ni_p);
1186: gp_x = 0.0;
1187: gp_y = 0.0;
1188: for (n = 0; n < NODES_PER_EL; n++) {
1189: gp_x = gp_x+Ni_p[n]*el_coords[2*n];
1190: gp_y = gp_y+Ni_p[n]*el_coords[2*n+1];
1191: }
1192: element_props[j][i].gp_coords[2*p] = gp_x;
1193: element_props[j][i].gp_coords[2*p+1] = gp_y;
1194: }
1195: }
1196: }
1198: /* define the coefficients */
1199: PetscOptionsGetBool(NULL,NULL,"-use_gp_coords",&use_gp_coords,NULL);
1201: for (j = sj; j < sj+ny; j++) {
1202: for (i = si; i < si+nx; i++) {
1203: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1204: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1205: PetscReal coord_x,coord_y;
1207: if (coefficient_structure == 0) {
1208: PetscReal opts_eta0,opts_eta1,opts_xc;
1209: PetscInt opts_nz;
1211: opts_eta0 = 1.0;
1212: opts_eta1 = 1.0;
1213: opts_xc = 0.5;
1214: opts_nz = 1;
1216: PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1217: PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1218: PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1219: PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);
1221: for (p = 0; p < GAUSS_POINTS; p++) {
1222: coord_x = centroid_x;
1223: coord_y = centroid_y;
1224: if (use_gp_coords) {
1225: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1226: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1227: }
1229: element_props[j][i].eta[p] = opts_eta0;
1230: if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;
1232: element_props[j][i].fx[p] = 0.0;
1233: element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1234: }
1235: } else if (coefficient_structure == 1) { /* square sinker */
1236: PetscReal opts_eta0,opts_eta1,opts_dx,opts_dy;
1238: opts_eta0 = 1.0;
1239: opts_eta1 = 1.0;
1240: opts_dx = 0.50;
1241: opts_dy = 0.50;
1243: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1244: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1245: PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1246: PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
1248: for (p = 0; p < GAUSS_POINTS; p++) {
1249: coord_x = centroid_x;
1250: coord_y = centroid_y;
1251: if (use_gp_coords) {
1252: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1253: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1254: }
1256: element_props[j][i].eta[p] = opts_eta0;
1257: element_props[j][i].fx[p] = 0.0;
1258: element_props[j][i].fy[p] = 0.0;
1260: if ((coord_x > -0.5*opts_dx+0.5) && (coord_x < 0.5*opts_dx+0.5)) {
1261: if ((coord_y > -0.5*opts_dy+0.5) && (coord_y < 0.5*opts_dy+0.5)) {
1262: element_props[j][i].eta[p] = opts_eta1;
1263: element_props[j][i].fx[p] = 0.0;
1264: element_props[j][i].fy[p] = -1.0;
1265: }
1266: }
1267: }
1268: } else if (coefficient_structure == 2) { /* circular sinker */
1269: PetscReal opts_eta0,opts_eta1,opts_r,radius2;
1271: opts_eta0 = 1.0;
1272: opts_eta1 = 1.0;
1273: opts_r = 0.25;
1275: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1276: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1277: PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);
1279: for (p = 0; p < GAUSS_POINTS; p++) {
1280: coord_x = centroid_x;
1281: coord_y = centroid_y;
1282: if (use_gp_coords) {
1283: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1284: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1285: }
1287: element_props[j][i].eta[p] = opts_eta0;
1288: element_props[j][i].fx[p] = 0.0;
1289: element_props[j][i].fy[p] = 0.0;
1291: radius2 = (coord_x-0.5)*(coord_x-0.5)+(coord_y-0.5)*(coord_y-0.5);
1292: if (radius2 < opts_r*opts_r) {
1293: element_props[j][i].eta[p] = opts_eta1;
1294: element_props[j][i].fx[p] = 0.0;
1295: element_props[j][i].fy[p] = -1.0;
1296: }
1297: }
1298: } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1299: PetscReal opts_eta0,opts_eta1,opts_r,opts_dx,opts_dy,opts_c0x,opts_c0y,opts_s0x,opts_s0y,opts_phi,radius2;
1301: opts_eta0 = 1.0;
1302: opts_eta1 = 1.0;
1303: opts_r = 0.25;
1304: opts_c0x = 0.35; /* circle center */
1305: opts_c0y = 0.35;
1306: opts_s0x = 0.7; /* square center */
1307: opts_s0y = 0.7;
1308: opts_dx = 0.25;
1309: opts_dy = 0.25;
1310: opts_phi = 25;
1312: PetscOptionsGetReal(NULL,NULL,"-sinker_eta0",&opts_eta0,NULL);
1313: PetscOptionsGetReal(NULL,NULL,"-sinker_eta1",&opts_eta1,NULL);
1314: PetscOptionsGetReal(NULL,NULL,"-sinker_r",&opts_r,NULL);
1315: PetscOptionsGetReal(NULL,NULL,"-sinker_c0x",&opts_c0x,NULL);
1316: PetscOptionsGetReal(NULL,NULL,"-sinker_c0y",&opts_c0y,NULL);
1317: PetscOptionsGetReal(NULL,NULL,"-sinker_s0x",&opts_s0x,NULL);
1318: PetscOptionsGetReal(NULL,NULL,"-sinker_s0y",&opts_s0y,NULL);
1319: PetscOptionsGetReal(NULL,NULL,"-sinker_dx",&opts_dx,NULL);
1320: PetscOptionsGetReal(NULL,NULL,"-sinker_dy",&opts_dy,NULL);
1321: PetscOptionsGetReal(NULL,NULL,"-sinker_phi",&opts_phi,NULL);
1322: opts_phi *= PETSC_PI / 180;
1324: for (p = 0; p < GAUSS_POINTS; p++) {
1325: coord_x = centroid_x;
1326: coord_y = centroid_y;
1327: if (use_gp_coords) {
1328: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1329: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1330: }
1332: element_props[j][i].eta[p] = opts_eta0;
1333: element_props[j][i].fx[p] = 0.0;
1334: element_props[j][i].fy[p] = -0.2;
1336: radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1337: if (radius2 < opts_r*opts_r
1338: || (PetscAbs(+(coord_x - opts_s0x)*PetscCosReal(opts_phi) + (coord_y - opts_s0y)*PetscSinReal(opts_phi)) < opts_dx/2
1339: && PetscAbs(-(coord_x - opts_s0x)*PetscSinReal(opts_phi) + (coord_y - opts_s0y)*PetscCosReal(opts_phi)) < opts_dy/2)) {
1340: element_props[j][i].eta[p] = opts_eta1;
1341: element_props[j][i].fx[p] = 0.0;
1342: element_props[j][i].fy[p] = -1.0;
1343: }
1344: }
1345: } else if (coefficient_structure == 4) { /* subdomain jump */
1346: PetscReal opts_mag,opts_eta0;
1347: PetscInt opts_nz,px,py;
1348: PetscBool jump;
1350: opts_mag = 1.0;
1351: opts_eta0 = 1.0;
1352: opts_nz = 1;
1354: PetscOptionsGetReal(NULL,NULL,"-jump_eta0",&opts_eta0,NULL);
1355: PetscOptionsGetReal(NULL,NULL,"-jump_magnitude",&opts_mag,NULL);
1356: PetscOptionsGetInt(NULL,NULL,"-jump_nz",&opts_nz,NULL);
1357: DMDAGetInfo(da_Stokes,NULL,NULL,NULL,NULL,&px,&py,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
1358: if (px%2) {
1359: jump = (PetscBool)(PetscGlobalRank%2);
1360: } else {
1361: jump = (PetscBool)((PetscGlobalRank/px)%2 ? PetscGlobalRank%2 : !(PetscGlobalRank%2));
1362: }
1363: for (p = 0; p < GAUSS_POINTS; p++) {
1364: coord_x = centroid_x;
1365: coord_y = centroid_y;
1366: if (use_gp_coords) {
1367: coord_x = PetscRealPart(element_props[j][i].gp_coords[2*p]);
1368: coord_y = PetscRealPart(element_props[j][i].gp_coords[2*p+1]);
1369: }
1371: element_props[j][i].eta[p] = jump ? PetscPowReal(10.0,opts_mag) : opts_eta0;
1372: element_props[j][i].fx[p] = 0.0;
1373: element_props[j][i].fy[p] = PetscSinReal(opts_nz*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1374: }
1375: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Unknown coefficient_structure");
1376: }
1377: }
1378: DMDAVecRestoreArray(prop_cda,prop_coords,&_prop_coords);
1380: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1382: DMDAVecRestoreArray(da_prop,l_properties,&element_props);
1383: DMLocalToGlobalBegin(da_prop,l_properties,ADD_VALUES,properties);
1384: DMLocalToGlobalEnd(da_prop,l_properties,ADD_VALUES,properties);
1386: if (output_gnuplot) {
1387: DMDACoordViewGnuplot2d(da_Stokes,"mesh");
1388: DMDAViewCoefficientsGnuplot2d(da_prop,properties,"Coefficients for Stokes eqn.","properties");
1389: }
1391: if (glvis) {
1392: Vec glv_prop,etaf;
1393: PetscViewer view;
1394: PetscInt dim;
1395: const char *fec = {"FiniteElementCollection: L2_2D_P1"};
1397: DMGetDimension(da_Stokes,&dim);
1398: VecCreateSeq(PETSC_COMM_SELF,GAUSS_POINTS*mx*mx,&etaf);
1399: PetscObjectSetName((PetscObject)etaf,"viscosity");
1400: PetscViewerGLVisOpen(PETSC_COMM_WORLD,PETSC_VIEWER_GLVIS_SOCKET,NULL,PETSC_DECIDE,&view);
1401: PetscViewerGLVisSetFields(view,1,&fec,&dim,glvis_extract_eta,(PetscObject*)&etaf,da_prop,NULL);
1402: DMGetLocalVector(da_prop,&glv_prop);
1403: DMGlobalToLocalBegin(da_prop,properties,INSERT_VALUES,glv_prop);
1404: DMGlobalToLocalEnd(da_prop,properties,INSERT_VALUES,glv_prop);
1405: VecSetDM(etaf,da_Stokes);
1406: VecView(glv_prop,view);
1407: PetscViewerDestroy(&view);
1408: DMRestoreLocalVector(da_prop,&glv_prop);
1409: VecDestroy(&etaf);
1410: }
1412: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1413: DMCreateMatrix(da_Stokes,&A);
1414: DMCreateMatrix(da_Stokes,&B);
1415: DMCreateGlobalVector(da_Stokes,&f);
1416: DMCreateGlobalVector(da_Stokes,&X);
1418: /* assemble A11 */
1419: MatZeroEntries(A);
1420: MatZeroEntries(B);
1421: VecZeroEntries(f);
1423: AssembleA_Stokes(A,da_Stokes,da_prop,properties);
1424: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
1425: /* build force vector */
1426: AssembleF_Stokes(f,da_Stokes,da_prop,properties);
1428: DMDABCApplyFreeSlip(da_Stokes,A,f);
1429: DMDABCApplyFreeSlip(da_Stokes,B,NULL);
1431: /* SOLVE */
1432: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1433: KSPSetOptionsPrefix(ksp_S,"stokes_");
1434: KSPSetDM(ksp_S,da_Stokes);
1435: KSPSetDMActive(ksp_S,PETSC_FALSE);
1436: KSPSetOperators(ksp_S,A,B);
1437: KSPSetFromOptions(ksp_S);
1438: {
1439: PC pc;
1440: const PetscInt ufields[] = {0,1},pfields[1] = {2};
1442: KSPGetPC(ksp_S,&pc);
1443: PCFieldSplitSetBlockSize(pc,3);
1444: PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1445: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1446: }
1448: {
1449: PC pc;
1450: PetscBool same = PETSC_FALSE;
1451: KSPGetPC(ksp_S,&pc);
1452: PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&same);
1453: if (same) {
1454: PetscBool usedivmat = PETSC_FALSE;
1455: KSPSetOperators(ksp_S,A,A);
1457: PetscOptionsGetBool(NULL,NULL,"-stokes_pc_bddc_use_divergence",&usedivmat,NULL);
1458: if (usedivmat) {
1459: IS *fields,vel;
1460: PetscInt i,nf;
1462: DMCreateFieldDecomposition(da_Stokes,&nf,NULL,&fields,NULL);
1463: ISConcatenate(PETSC_COMM_WORLD,2,fields,&vel);
1464: MatZeroRowsIS(B,fields[2],1.0,NULL,NULL); /* we put 1.0 on the diagonal to pick the pressure average too */
1465: MatTranspose(B,MAT_INPLACE_MATRIX,&B);
1466: MatZeroRowsIS(B,vel,0.0,NULL,NULL);
1467: ISDestroy(&vel);
1468: PCBDDCSetDivergenceMat(pc,B,PETSC_FALSE,NULL);
1469: for (i=0;i<nf;i++) {
1470: ISDestroy(&fields[i]);
1471: }
1472: PetscFree(fields);
1473: }
1474: }
1475: }
1477: {
1478: PC pc_S;
1479: KSP *sub_ksp,ksp_U;
1480: PetscInt nsplits;
1481: DM da_U;
1482: PetscBool is_pcfs;
1484: KSPSetUp(ksp_S);
1485: KSPGetPC(ksp_S,&pc_S);
1487: is_pcfs = PETSC_FALSE;
1488: PetscObjectTypeCompare((PetscObject)pc_S,PCFIELDSPLIT,&is_pcfs);
1490: if (is_pcfs) {
1491: PCFieldSplitGetSubKSP(pc_S,&nsplits,&sub_ksp);
1492: ksp_U = sub_ksp[0];
1493: PetscFree(sub_ksp);
1495: if (nsplits == 2) {
1496: DMDACreateCompatibleDMDA(da_Stokes,2,&da_U);
1498: KSPSetDM(ksp_U,da_U);
1499: KSPSetDMActive(ksp_U,PETSC_FALSE);
1500: DMDestroy(&da_U);
1501: }
1502: }
1503: }
1505: KSPSolve(ksp_S,f,X);
1507: PetscOptionsGetString(NULL,NULL,"-o",filename,sizeof(filename),&set);
1508: if (set) {
1509: char *ext;
1510: PetscViewer viewer;
1511: PetscBool flg;
1512: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1513: PetscStrrchr(filename,'.',&ext);
1514: PetscStrcmp("vts",ext,&flg);
1515: if (flg) {
1516: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1517: } else {
1518: PetscViewerSetType(viewer,PETSCVIEWERBINARY);
1519: }
1520: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1521: PetscViewerFileSetName(viewer,filename);
1522: VecView(X,viewer);
1523: PetscViewerDestroy(&viewer);
1524: }
1525: if (output_gnuplot) {
1526: DMDAViewGnuplot2d(da_Stokes,X,"Velocity solution for Stokes eqn.","X");
1527: }
1529: if (glvis) {
1530: PetscViewer view;
1532: PetscViewerCreate(PETSC_COMM_WORLD,&view);
1533: PetscViewerSetType(view,PETSCVIEWERGLVIS);
1534: VecView(X,view);
1535: PetscViewerDestroy(&view);
1536: }
1538: KSPGetIterationNumber(ksp_S,&its);
1540: if (coefficient_structure == 0) {
1541: PetscReal opts_eta0,opts_eta1,opts_xc;
1542: PetscInt opts_nz,N;
1543: DM da_Stokes_analytic;
1544: Vec X_analytic;
1545: PetscReal nrm1[3],nrm2[3],nrmI[3];
1547: opts_eta0 = 1.0;
1548: opts_eta1 = 1.0;
1549: opts_xc = 0.5;
1550: opts_nz = 1;
1552: PetscOptionsGetReal(NULL,NULL,"-solcx_eta0",&opts_eta0,NULL);
1553: PetscOptionsGetReal(NULL,NULL,"-solcx_eta1",&opts_eta1,NULL);
1554: PetscOptionsGetReal(NULL,NULL,"-solcx_xc",&opts_xc,NULL);
1555: PetscOptionsGetInt(NULL,NULL,"-solcx_nz",&opts_nz,NULL);
1557: DMDACreateSolCx(opts_eta0,opts_eta1,opts_xc,opts_nz,mx,my,&da_Stokes_analytic,&X_analytic);
1558: if (output_gnuplot) {
1559: DMDAViewGnuplot2d(da_Stokes_analytic,X_analytic,"Analytic solution for Stokes eqn.","X_analytic");
1560: }
1561: DMDAIntegrateErrors(da_Stokes_analytic,X,X_analytic);
1563: VecAXPY(X_analytic,-1.0,X);
1564: VecGetSize(X_analytic,&N);
1565: N = N/3;
1567: VecStrideNorm(X_analytic,0,NORM_1,&nrm1[0]);
1568: VecStrideNorm(X_analytic,0,NORM_2,&nrm2[0]);
1569: VecStrideNorm(X_analytic,0,NORM_INFINITY,&nrmI[0]);
1571: VecStrideNorm(X_analytic,1,NORM_1,&nrm1[1]);
1572: VecStrideNorm(X_analytic,1,NORM_2,&nrm2[1]);
1573: VecStrideNorm(X_analytic,1,NORM_INFINITY,&nrmI[1]);
1575: VecStrideNorm(X_analytic,2,NORM_1,&nrm1[2]);
1576: VecStrideNorm(X_analytic,2,NORM_2,&nrm2[2]);
1577: VecStrideNorm(X_analytic,2,NORM_INFINITY,&nrmI[2]);
1579: DMDestroy(&da_Stokes_analytic);
1580: VecDestroy(&X_analytic);
1581: }
1583: KSPDestroy(&ksp_S);
1584: VecDestroy(&X);
1585: VecDestroy(&f);
1586: MatDestroy(&A);
1587: MatDestroy(&B);
1589: DMDestroy(&da_Stokes);
1590: DMDestroy(&da_prop);
1592: VecDestroy(&properties);
1593: VecDestroy(&l_properties);
1594: return(0);
1595: }
1597: int main(int argc,char **args)
1598: {
1600: PetscInt mx,my;
1602: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
1603: mx = my = 10;
1604: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1605: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1606: solve_stokes_2d_coupled(mx,my);
1607: PetscFinalize();
1608: return ierr;
1609: }
1611: /* -------------------------- helpers for boundary conditions -------------------------------- */
1612: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
1613: {
1614: DM cda;
1615: Vec coords;
1616: PetscInt si,sj,nx,ny,i,j;
1617: PetscInt M,N;
1618: DMDACoor2d **_coords;
1619: const PetscInt *g_idx;
1620: PetscInt *bc_global_ids;
1621: PetscScalar *bc_vals;
1622: PetscInt nbcs;
1623: PetscInt n_dofs;
1624: PetscErrorCode ierr;
1625: ISLocalToGlobalMapping ltogm;
1628: DMGetLocalToGlobalMapping(da,<ogm);
1629: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1631: DMGetCoordinateDM(da,&cda);
1632: DMGetCoordinatesLocal(da,&coords);
1633: DMDAVecGetArray(cda,coords,&_coords);
1634: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1635: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1637: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1638: PetscMalloc1(ny*n_dofs,&bc_vals);
1640: /* init the entries to -1 so VecSetValues will ignore them */
1641: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1643: i = nx-1;
1644: for (j = 0; j < ny; j++) {
1645: PetscInt local_id;
1647: local_id = i+j*nx;
1649: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1651: bc_vals[j] = 0.0;
1652: }
1653: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1654: nbcs = 0;
1655: if ((si+nx) == (M)) nbcs = ny;
1657: if (b) {
1658: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1659: VecAssemblyBegin(b);
1660: VecAssemblyEnd(b);
1661: }
1662: if (A) {
1663: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1664: }
1666: PetscFree(bc_vals);
1667: PetscFree(bc_global_ids);
1669: DMDAVecRestoreArray(cda,coords,&_coords);
1670: return(0);
1671: }
1673: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
1674: {
1675: DM cda;
1676: Vec coords;
1677: PetscInt si,sj,nx,ny,i,j;
1678: PetscInt M,N;
1679: DMDACoor2d **_coords;
1680: const PetscInt *g_idx;
1681: PetscInt *bc_global_ids;
1682: PetscScalar *bc_vals;
1683: PetscInt nbcs;
1684: PetscInt n_dofs;
1685: PetscErrorCode ierr;
1686: ISLocalToGlobalMapping ltogm;
1689: DMGetLocalToGlobalMapping(da,<ogm);
1690: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1692: DMGetCoordinateDM(da,&cda);
1693: DMGetCoordinatesLocal(da,&coords);
1694: DMDAVecGetArray(cda,coords,&_coords);
1695: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1696: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1698: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1699: PetscMalloc1(ny*n_dofs,&bc_vals);
1701: /* init the entries to -1 so VecSetValues will ignore them */
1702: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1704: i = 0;
1705: for (j = 0; j < ny; j++) {
1706: PetscInt local_id;
1708: local_id = i+j*nx;
1710: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1712: bc_vals[j] = 0.0;
1713: }
1714: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1715: nbcs = 0;
1716: if (si == 0) nbcs = ny;
1718: if (b) {
1719: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1720: VecAssemblyBegin(b);
1721: VecAssemblyEnd(b);
1722: }
1724: if (A) {
1725: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1726: }
1728: PetscFree(bc_vals);
1729: PetscFree(bc_global_ids);
1731: DMDAVecRestoreArray(cda,coords,&_coords);
1732: return(0);
1733: }
1735: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
1736: {
1737: DM cda;
1738: Vec coords;
1739: PetscInt si,sj,nx,ny,i,j;
1740: PetscInt M,N;
1741: DMDACoor2d **_coords;
1742: const PetscInt *g_idx;
1743: PetscInt *bc_global_ids;
1744: PetscScalar *bc_vals;
1745: PetscInt nbcs;
1746: PetscInt n_dofs;
1747: PetscErrorCode ierr;
1748: ISLocalToGlobalMapping ltogm;
1751: DMGetLocalToGlobalMapping(da,<ogm);
1752: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1754: DMGetCoordinateDM(da,&cda);
1755: DMGetCoordinatesLocal(da,&coords);
1756: DMDAVecGetArray(cda,coords,&_coords);
1757: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1758: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1760: PetscMalloc1(nx,&bc_global_ids);
1761: PetscMalloc1(nx,&bc_vals);
1763: /* init the entries to -1 so VecSetValues will ignore them */
1764: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1766: j = ny-1;
1767: for (i = 0; i < nx; i++) {
1768: PetscInt local_id;
1770: local_id = i+j*nx;
1772: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1774: bc_vals[i] = 0.0;
1775: }
1776: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1777: nbcs = 0;
1778: if ((sj+ny) == (N)) nbcs = nx;
1780: if (b) {
1781: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1782: VecAssemblyBegin(b);
1783: VecAssemblyEnd(b);
1784: }
1785: if (A) {
1786: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1787: }
1789: PetscFree(bc_vals);
1790: PetscFree(bc_global_ids);
1792: DMDAVecRestoreArray(cda,coords,&_coords);
1793: return(0);
1794: }
1796: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
1797: {
1798: DM cda;
1799: Vec coords;
1800: PetscInt si,sj,nx,ny,i,j;
1801: PetscInt M,N;
1802: DMDACoor2d **_coords;
1803: const PetscInt *g_idx;
1804: PetscInt *bc_global_ids;
1805: PetscScalar *bc_vals;
1806: PetscInt nbcs;
1807: PetscInt n_dofs;
1808: PetscErrorCode ierr;
1809: ISLocalToGlobalMapping ltogm;
1812: DMGetLocalToGlobalMapping(da,<ogm);
1813: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1815: DMGetCoordinateDM(da,&cda);
1816: DMGetCoordinatesLocal(da,&coords);
1817: DMDAVecGetArray(cda,coords,&_coords);
1818: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1819: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1821: PetscMalloc1(nx,&bc_global_ids);
1822: PetscMalloc1(nx,&bc_vals);
1824: /* init the entries to -1 so VecSetValues will ignore them */
1825: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1827: j = 0;
1828: for (i = 0; i < nx; i++) {
1829: PetscInt local_id;
1831: local_id = i+j*nx;
1833: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1835: bc_vals[i] = 0.0;
1836: }
1837: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1838: nbcs = 0;
1839: if (sj == 0) nbcs = nx;
1841: if (b) {
1842: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1843: VecAssemblyBegin(b);
1844: VecAssemblyEnd(b);
1845: }
1846: if (A) {
1847: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1848: }
1850: PetscFree(bc_vals);
1851: PetscFree(bc_global_ids);
1853: DMDAVecRestoreArray(cda,coords,&_coords);
1854: return(0);
1855: }
1857: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
1858: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)
1859: {
1863: BCApplyZero_NORTH(da_Stokes,1,A,f);
1864: BCApplyZero_EAST(da_Stokes,0,A,f);
1865: BCApplyZero_SOUTH(da_Stokes,1,A,f);
1866: BCApplyZero_WEST(da_Stokes,0,A,f);
1867: return(0);
1868: }
1871: /*TEST
1873: build:
1874: requires: !complex !single
1876: test:
1877: 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
1879: test:
1880: suffix: 2
1881: 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
1882: output_file: output/ex43_1.out
1884: test:
1885: suffix: 3
1886: nsize: 4
1887: 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
1889: test:
1890: suffix: 4
1891: nsize: 4
1892: 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
1894: test:
1895: suffix: 5
1896: nsize: 4
1897: 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
1899: test:
1900: suffix: 6
1901: nsize: 8
1902: 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
1904: test:
1905: suffix: bjacobi
1906: nsize: 4
1907: args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type aij -stokes_ksp_converged_reason
1909: test:
1910: suffix: bjacobi_baij
1911: nsize: 4
1912: args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type baij -stokes_ksp_converged_reason
1913: output_file: output/ex43_bjacobi.out
1915: test:
1916: suffix: nested_gmg
1917: nsize: 4
1918: 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
1920: test:
1921: suffix: fetidp
1922: nsize: 8
1923: 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
1925: test:
1926: suffix: fetidp_unsym
1927: nsize: 8
1928: 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
1930: test:
1931: suffix: bddc_stokes_deluxe
1932: nsize: 8
1933: 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
1935: test:
1936: suffix: bddc_stokes_subdomainjump_deluxe
1937: nsize: 9
1938: 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
1941: TEST*/