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 evaluated 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++) array[cum++] = props[ej][ei].eta[p];
102: }
103: }
104: VecRestoreArray(Vf[0], &array);
105: DMDAVecRestoreArrayRead(properties_da, V, &props);
106: return 0;
107: }
109: /*
110: Element: Local basis function ordering
111: 1-----2
112: | |
113: | |
114: 0-----3
115: */
116: static void ConstructQ12D_Ni(PetscScalar _xi[], PetscScalar Ni[])
117: {
118: PetscScalar xi = _xi[0];
119: PetscScalar eta = _xi[1];
121: Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
122: Ni[1] = 0.25 * (1.0 - xi) * (1.0 + eta);
123: Ni[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
124: Ni[3] = 0.25 * (1.0 + xi) * (1.0 - eta);
125: }
127: static void ConstructQ12D_GNi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
128: {
129: PetscScalar xi = _xi[0];
130: PetscScalar eta = _xi[1];
132: GNi[0][0] = -0.25 * (1.0 - eta);
133: GNi[0][1] = -0.25 * (1.0 + eta);
134: GNi[0][2] = 0.25 * (1.0 + eta);
135: GNi[0][3] = 0.25 * (1.0 - eta);
137: GNi[1][0] = -0.25 * (1.0 - xi);
138: GNi[1][1] = 0.25 * (1.0 - xi);
139: GNi[1][2] = 0.25 * (1.0 + xi);
140: GNi[1][3] = -0.25 * (1.0 + xi);
141: }
143: static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL], PetscScalar GNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
144: {
145: PetscScalar J00, J01, J10, J11, J;
146: PetscScalar iJ00, iJ01, iJ10, iJ11;
147: PetscInt i;
149: J00 = J01 = J10 = J11 = 0.0;
150: for (i = 0; i < NODES_PER_EL; i++) {
151: PetscScalar cx = coords[2 * i];
152: PetscScalar cy = coords[2 * i + 1];
154: J00 = J00 + GNi[0][i] * cx; /* J_xx = dx/dxi */
155: J01 = J01 + GNi[0][i] * cy; /* J_xy = dy/dxi */
156: J10 = J10 + GNi[1][i] * cx; /* J_yx = dx/deta */
157: J11 = J11 + GNi[1][i] * cy; /* J_yy = dy/deta */
158: }
159: J = (J00 * J11) - (J01 * J10);
161: iJ00 = J11 / J;
162: iJ01 = -J01 / J;
163: iJ10 = -J10 / J;
164: iJ11 = J00 / J;
166: for (i = 0; i < NODES_PER_EL; i++) {
167: GNx[0][i] = GNi[0][i] * iJ00 + GNi[1][i] * iJ01;
168: GNx[1][i] = GNi[0][i] * iJ10 + GNi[1][i] * iJ11;
169: }
171: if (det_J) *det_J = J;
172: }
174: static void ConstructGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
175: {
176: *ngp = 4;
177: gp_xi[0][0] = -0.57735026919;
178: gp_xi[0][1] = -0.57735026919;
179: gp_xi[1][0] = -0.57735026919;
180: gp_xi[1][1] = 0.57735026919;
181: gp_xi[2][0] = 0.57735026919;
182: gp_xi[2][1] = 0.57735026919;
183: gp_xi[3][0] = 0.57735026919;
184: 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;
202: s_u[0].j = j;
203: s_u[0].c = 0; /* Vx0 */
204: s_u[1].i = i;
205: s_u[1].j = j;
206: s_u[1].c = 1; /* Vy0 */
208: /* node 1 */
209: s_u[2].i = i;
210: s_u[2].j = j + 1;
211: s_u[2].c = 0; /* Vx1 */
212: s_u[3].i = i;
213: s_u[3].j = j + 1;
214: s_u[3].c = 1; /* Vy1 */
216: /* node 2 */
217: s_u[4].i = i + 1;
218: s_u[4].j = j + 1;
219: s_u[4].c = 0; /* Vx2 */
220: s_u[5].i = i + 1;
221: s_u[5].j = j + 1;
222: s_u[5].c = 1; /* Vy2 */
224: /* node 3 */
225: s_u[6].i = i + 1;
226: s_u[6].j = j;
227: s_u[6].c = 0; /* Vx3 */
228: s_u[7].i = i + 1;
229: s_u[7].j = j;
230: s_u[7].c = 1; /* Vy3 */
232: /* pressure */
233: s_p[0].i = i;
234: s_p[0].j = j;
235: s_p[0].c = 2; /* P0 */
236: s_p[1].i = i;
237: s_p[1].j = j + 1;
238: s_p[1].c = 2; /* P1 */
239: s_p[2].i = i + 1;
240: s_p[2].j = j + 1;
241: s_p[2].c = 2; /* P2 */
242: s_p[3].i = i + 1;
243: s_p[3].j = j;
244: s_p[3].c = 2; /* P3 */
245: return 0;
246: }
248: static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da, PetscInt **_lx, PetscInt **_ly)
249: {
250: PetscMPIInt rank;
251: PetscInt proc_I, proc_J;
252: PetscInt cpu_x, cpu_y;
253: PetscInt local_mx, local_my;
254: Vec vlx, vly;
255: PetscInt *LX, *LY, i;
256: PetscScalar *_a;
257: Vec V_SEQ;
258: VecScatter ctx;
261: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
263: DMDAGetInfo(da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0);
265: proc_J = rank / cpu_x;
266: proc_I = rank - cpu_x * proc_J;
268: PetscMalloc1(cpu_x, &LX);
269: PetscMalloc1(cpu_y, &LY);
271: DMDAGetElementsSizes(da, &local_mx, &local_my, NULL);
272: VecCreate(PETSC_COMM_WORLD, &vlx);
273: VecSetSizes(vlx, PETSC_DECIDE, cpu_x);
274: VecSetFromOptions(vlx);
276: VecCreate(PETSC_COMM_WORLD, &vly);
277: VecSetSizes(vly, PETSC_DECIDE, cpu_y);
278: VecSetFromOptions(vly);
280: VecSetValue(vlx, proc_I, (PetscScalar)(local_mx + 1.0e-9), INSERT_VALUES);
281: VecSetValue(vly, proc_J, (PetscScalar)(local_my + 1.0e-9), INSERT_VALUES);
282: VecAssemblyBegin(vlx); VecAssemblyEnd(vlx);
283: VecAssemblyBegin(vly); VecAssemblyEnd(vly);
285: VecScatterCreateToAll(vlx, &ctx, &V_SEQ);
286: VecScatterBegin(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD);
287: VecScatterEnd(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD);
288: VecGetArray(V_SEQ, &_a);
289: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
290: VecRestoreArray(V_SEQ, &_a);
291: VecScatterDestroy(&ctx);
292: VecDestroy(&V_SEQ);
294: VecScatterCreateToAll(vly, &ctx, &V_SEQ);
295: VecScatterBegin(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD);
296: VecScatterEnd(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD);
297: VecGetArray(V_SEQ, &_a);
298: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
299: VecRestoreArray(V_SEQ, &_a);
300: VecScatterDestroy(&ctx);
301: VecDestroy(&V_SEQ);
303: *_lx = LX;
304: *_ly = LY;
306: VecDestroy(&vlx);
307: VecDestroy(&vly);
308: return 0;
309: }
311: static PetscErrorCode DMDACoordViewGnuplot2d(DM da, const char prefix[])
312: {
313: DM cda;
314: Vec coords;
315: DMDACoor2d **_coords;
316: PetscInt si, sj, nx, ny, i, j;
317: FILE *fp;
318: char fname[PETSC_MAX_PATH_LEN];
319: PetscMPIInt rank;
322: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
323: PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank);
324: PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp);
327: PetscFPrintf(PETSC_COMM_SELF, fp, "### Element geometry for processor %1.4d ### \n", rank);
329: DMGetCoordinateDM(da, &cda);
330: DMGetCoordinatesLocal(da, &coords);
331: DMDAVecGetArrayRead(cda, coords, &_coords);
332: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
333: for (j = sj; j < sj + ny - 1; j++) {
334: for (i = si; i < si + nx - 1; i++) {
335: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y));
336: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i].x), (double)PetscRealPart(_coords[j + 1][i].y));
337: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i + 1].x), (double)PetscRealPart(_coords[j + 1][i + 1].y));
338: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i + 1].x), (double)PetscRealPart(_coords[j][i + 1].y));
339: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n\n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y));
340: }
341: }
342: DMDAVecRestoreArrayRead(cda, coords, &_coords);
344: PetscFClose(PETSC_COMM_SELF, fp);
345: return 0;
346: }
348: static PetscErrorCode DMDAViewGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
349: {
350: DM cda;
351: Vec coords, local_fields;
352: DMDACoor2d **_coords;
353: FILE *fp;
354: char fname[PETSC_MAX_PATH_LEN];
355: PetscMPIInt rank;
356: PetscInt si, sj, nx, ny, i, j;
357: PetscInt n_dofs, d;
358: PetscScalar *_fields;
361: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
362: PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank);
363: PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp);
366: PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank);
367: DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);
368: PetscFPrintf(PETSC_COMM_SELF, fp, "### x y ");
369: for (d = 0; d < n_dofs; d++) {
370: const char *field_name;
371: DMDAGetFieldName(da, d, &field_name);
372: PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name);
373: }
374: PetscFPrintf(PETSC_COMM_SELF, fp, "###\n");
376: DMGetCoordinateDM(da, &cda);
377: DMGetCoordinatesLocal(da, &coords);
378: DMDAVecGetArray(cda, coords, &_coords);
379: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
381: DMCreateLocalVector(da, &local_fields);
382: DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields);
383: DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields);
384: VecGetArray(local_fields, &_fields);
386: for (j = sj; j < sj + ny; j++) {
387: for (i = si; i < si + nx; i++) {
388: PetscScalar coord_x, coord_y;
389: PetscScalar field_d;
391: coord_x = _coords[j][i].x;
392: coord_y = _coords[j][i].y;
394: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y));
395: for (d = 0; d < n_dofs; d++) {
396: field_d = _fields[n_dofs * ((i - si) + (j - sj) * (nx)) + d];
397: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e ", (double)PetscRealPart(field_d));
398: }
399: PetscFPrintf(PETSC_COMM_SELF, fp, "\n");
400: }
401: }
402: VecRestoreArray(local_fields, &_fields);
403: VecDestroy(&local_fields);
405: DMDAVecRestoreArray(cda, coords, &_coords);
407: PetscFClose(PETSC_COMM_SELF, fp);
408: return 0;
409: }
411: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
412: {
413: DM cda;
414: Vec local_fields;
415: FILE *fp;
416: char fname[PETSC_MAX_PATH_LEN];
417: PetscMPIInt rank;
418: PetscInt si, sj, nx, ny, i, j, p;
419: PetscInt n_dofs, d;
420: GaussPointCoefficients **_coefficients;
423: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
424: PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank);
425: PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp);
428: PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank);
429: DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);
430: PetscFPrintf(PETSC_COMM_SELF, fp, "### x y ");
431: for (d = 0; d < n_dofs; d++) {
432: const char *field_name;
433: DMDAGetFieldName(da, d, &field_name);
434: PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name);
435: }
436: PetscFPrintf(PETSC_COMM_SELF, fp, "###\n");
438: DMGetCoordinateDM(da, &cda);
439: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
441: DMCreateLocalVector(da, &local_fields);
442: DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields);
443: DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields);
444: DMDAVecGetArray(da, local_fields, &_coefficients);
446: for (j = sj; j < sj + ny; j++) {
447: for (i = si; i < si + nx; i++) {
448: PetscScalar coord_x, coord_y;
450: for (p = 0; p < GAUSS_POINTS; p++) {
451: coord_x = _coefficients[j][i].gp_coords[2 * p];
452: coord_y = _coefficients[j][i].gp_coords[2 * p + 1];
454: PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y));
456: 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]));
457: PetscFPrintf(PETSC_COMM_SELF, fp, "\n");
458: }
459: }
460: }
461: DMDAVecRestoreArray(da, local_fields, &_coefficients);
462: VecDestroy(&local_fields);
464: PetscFClose(PETSC_COMM_SELF, fp);
465: return 0;
466: }
468: 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)
469: {
470: PetscInt ij;
471: PetscInt r, c, nc;
473: nc = u_NPE * u_dof;
474: r = w_dof * wi + wd;
475: c = u_dof * ui + ud;
476: ij = r * nc + c;
477: return ij;
478: }
480: /*
481: D = [ 2.eta 0 0 ]
482: [ 0 2.eta 0 ]
483: [ 0 0 eta ]
485: B = [ d_dx 0 ]
486: [ 0 d_dy ]
487: [ d_dy d_dx ]
488: */
489: static void FormStressOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
490: {
491: PetscInt ngp;
492: PetscScalar gp_xi[GAUSS_POINTS][2];
493: PetscScalar gp_weight[GAUSS_POINTS];
494: PetscInt p, i, j, k;
495: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
496: PetscScalar J_p, tildeD[3];
497: PetscScalar B[3][U_DOFS * NODES_PER_EL];
499: /* define quadrature rule */
500: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
502: /* evaluate integral */
503: for (p = 0; p < ngp; p++) {
504: ConstructQ12D_GNi(gp_xi[p], GNi_p);
505: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
507: for (i = 0; i < NODES_PER_EL; i++) {
508: PetscScalar d_dx_i = GNx_p[0][i];
509: PetscScalar d_dy_i = GNx_p[1][i];
511: B[0][2 * i] = d_dx_i;
512: B[0][2 * i + 1] = 0.0;
513: B[1][2 * i] = 0.0;
514: B[1][2 * i + 1] = d_dy_i;
515: B[2][2 * i] = d_dy_i;
516: B[2][2 * i + 1] = d_dx_i;
517: }
519: tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
520: tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
521: tildeD[2] = gp_weight[p] * J_p * eta[p];
523: /* form Bt tildeD B */
524: /*
525: Ke_ij = Bt_ik . D_kl . B_lj
526: = B_ki . D_kl . B_lj
527: = B_ki . D_kk . B_kj
528: */
529: for (i = 0; i < 8; i++) {
530: for (j = 0; j < 8; j++) {
531: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
532: Ke[i + 8 * j] = Ke[i + 8 * j] + B[k][i] * tildeD[k] * B[k][j];
533: }
534: }
535: }
536: }
537: }
539: static void FormGradientOperatorQ1(PetscScalar Ke[], PetscScalar coords[])
540: {
541: PetscInt ngp;
542: PetscScalar gp_xi[GAUSS_POINTS][2];
543: PetscScalar gp_weight[GAUSS_POINTS];
544: PetscInt p, i, j, di;
545: PetscScalar Ni_p[NODES_PER_EL];
546: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
547: PetscScalar J_p, fac;
549: /* define quadrature rule */
550: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
552: /* evaluate integral */
553: for (p = 0; p < ngp; p++) {
554: ConstructQ12D_Ni(gp_xi[p], Ni_p);
555: ConstructQ12D_GNi(gp_xi[p], GNi_p);
556: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
557: fac = gp_weight[p] * J_p;
559: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
560: for (di = 0; di < NSD; di++) { /* u dofs */
561: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
562: PetscInt IJ;
563: IJ = ASS_MAP_wIwDI_uJuDJ(i, di, NODES_PER_EL, 2, j, 0, NODES_PER_EL, 1);
565: Ke[IJ] = Ke[IJ] - GNx_p[di][i] * Ni_p[j] * fac;
566: }
567: }
568: }
569: }
570: }
572: static void FormDivergenceOperatorQ1(PetscScalar De[], PetscScalar coords[])
573: {
574: PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
575: PetscInt i, j;
576: PetscInt nr_g, nc_g;
578: PetscMemzero(Ge, sizeof(Ge));
579: FormGradientOperatorQ1(Ge, coords);
581: nr_g = U_DOFS * NODES_PER_EL;
582: nc_g = P_DOFS * NODES_PER_EL;
584: for (i = 0; i < nr_g; i++) {
585: for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
586: }
587: }
589: static void FormStabilisationOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
590: {
591: PetscInt ngp;
592: PetscScalar gp_xi[GAUSS_POINTS][2];
593: PetscScalar gp_weight[GAUSS_POINTS];
594: PetscInt p, i, j;
595: PetscScalar Ni_p[NODES_PER_EL];
596: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
597: PetscScalar J_p, fac, eta_avg;
599: /* define quadrature rule */
600: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
602: /* evaluate integral */
603: for (p = 0; p < ngp; p++) {
604: ConstructQ12D_Ni(gp_xi[p], Ni_p);
605: ConstructQ12D_GNi(gp_xi[p], GNi_p);
606: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
607: fac = gp_weight[p] * J_p;
609: for (i = 0; i < NODES_PER_EL; i++) {
610: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = Ke[NODES_PER_EL * i + j] - fac * (Ni_p[i] * Ni_p[j] - 0.0625);
611: }
612: }
614: /* scale */
615: eta_avg = 0.0;
616: for (p = 0; p < ngp; p++) eta_avg += eta[p];
617: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
618: fac = 1.0 / eta_avg;
619: for (i = 0; i < NODES_PER_EL; i++) {
620: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
621: }
622: }
624: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
625: {
626: PetscInt ngp;
627: PetscScalar gp_xi[GAUSS_POINTS][2];
628: PetscScalar gp_weight[GAUSS_POINTS];
629: PetscInt p, i, j;
630: PetscScalar Ni_p[NODES_PER_EL];
631: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
632: PetscScalar J_p, fac, eta_avg;
634: /* define quadrature rule */
635: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
637: /* evaluate integral */
638: for (p = 0; p < ngp; p++) {
639: ConstructQ12D_Ni(gp_xi[p], Ni_p);
640: ConstructQ12D_GNi(gp_xi[p], GNi_p);
641: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
642: fac = gp_weight[p] * J_p;
644: for (i = 0; i < NODES_PER_EL; i++) {
645: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = Ke[NODES_PER_EL * i + j] - fac * Ni_p[i] * Ni_p[j];
646: }
647: }
649: /* scale */
650: eta_avg = 0.0;
651: for (p = 0; p < ngp; p++) eta_avg += eta[p];
652: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
653: fac = 1.0 / eta_avg;
654: for (i = 0; i < NODES_PER_EL; i++) {
655: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
656: }
657: }
659: static void FormMomentumRhsQ1(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
660: {
661: PetscInt ngp;
662: PetscScalar gp_xi[GAUSS_POINTS][2];
663: PetscScalar gp_weight[GAUSS_POINTS];
664: PetscInt p, i;
665: PetscScalar Ni_p[NODES_PER_EL];
666: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
667: PetscScalar J_p, fac;
669: /* define quadrature rule */
670: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
672: /* evaluate integral */
673: for (p = 0; p < ngp; p++) {
674: ConstructQ12D_Ni(gp_xi[p], Ni_p);
675: ConstructQ12D_GNi(gp_xi[p], GNi_p);
676: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
677: fac = gp_weight[p] * J_p;
679: for (i = 0; i < NODES_PER_EL; i++) {
680: Fe[NSD * i] += fac * Ni_p[i] * fx[p];
681: Fe[NSD * i + 1] += fac * Ni_p[i] * fy[p];
682: }
683: }
684: }
686: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords, PetscInt ei, PetscInt ej, PetscScalar el_coords[])
687: {
689: /* get coords for the element */
690: el_coords[NSD * 0] = _coords[ej][ei].x;
691: el_coords[NSD * 0 + 1] = _coords[ej][ei].y;
692: el_coords[NSD * 1] = _coords[ej + 1][ei].x;
693: el_coords[NSD * 1 + 1] = _coords[ej + 1][ei].y;
694: el_coords[NSD * 2] = _coords[ej + 1][ei + 1].x;
695: el_coords[NSD * 2 + 1] = _coords[ej + 1][ei + 1].y;
696: el_coords[NSD * 3] = _coords[ej][ei + 1].x;
697: el_coords[NSD * 3 + 1] = _coords[ej][ei + 1].y;
698: return 0;
699: }
701: static PetscErrorCode AssembleA_Stokes(Mat A, DM stokes_da, DM properties_da, Vec properties)
702: {
703: DM cda;
704: Vec coords;
705: DMDACoor2d **_coords;
706: MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
707: MatStencil p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
708: PetscInt sex, sey, mx, my;
709: PetscInt ei, ej;
710: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
711: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
712: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
713: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
714: PetscScalar el_coords[NODES_PER_EL * NSD];
715: Vec local_properties;
716: GaussPointCoefficients **props;
717: PetscScalar *prop_eta;
720: /* setup for coords */
721: DMGetCoordinateDM(stokes_da, &cda);
722: DMGetCoordinatesLocal(stokes_da, &coords);
723: DMDAVecGetArray(cda, coords, &_coords);
725: /* setup for coefficients */
726: DMCreateLocalVector(properties_da, &local_properties);
727: DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties);
728: DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties);
729: DMDAVecGetArray(properties_da, local_properties, &props);
731: DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL);
732: DMDAGetElementsSizes(stokes_da, &mx, &my, NULL);
733: for (ej = sey; ej < sey + my; ej++) {
734: for (ei = sex; ei < sex + mx; ei++) {
735: /* get coords for the element */
736: GetElementCoords(_coords, ei, ej, el_coords);
738: /* get coefficients for the element */
739: prop_eta = props[ej][ei].eta;
741: /* initialise element stiffness matrix */
742: PetscMemzero(Ae, sizeof(Ae));
743: PetscMemzero(Ge, sizeof(Ge));
744: PetscMemzero(De, sizeof(De));
745: PetscMemzero(Ce, sizeof(Ce));
747: /* form element stiffness matrix */
748: FormStressOperatorQ1(Ae, el_coords, prop_eta);
749: FormGradientOperatorQ1(Ge, el_coords);
750: FormDivergenceOperatorQ1(De, el_coords);
751: FormStabilisationOperatorQ1(Ce, el_coords, prop_eta);
753: /* insert element matrix into global matrix */
754: DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej);
755: MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES);
756: MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES);
757: MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES);
758: MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES);
759: }
760: }
761: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
762: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
764: DMDAVecRestoreArray(cda, coords, &_coords);
766: DMDAVecRestoreArray(properties_da, local_properties, &props);
767: VecDestroy(&local_properties);
768: return 0;
769: }
771: static PetscErrorCode AssembleA_PCStokes(Mat A, DM stokes_da, DM properties_da, Vec properties)
772: {
773: DM cda;
774: Vec coords;
775: DMDACoor2d **_coords;
776: MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
777: MatStencil p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
778: PetscInt sex, sey, mx, my;
779: PetscInt ei, ej;
780: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
781: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
782: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
783: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
784: PetscScalar el_coords[NODES_PER_EL * NSD];
785: Vec local_properties;
786: GaussPointCoefficients **props;
787: PetscScalar *prop_eta;
790: /* setup for coords */
791: DMGetCoordinateDM(stokes_da, &cda);
792: DMGetCoordinatesLocal(stokes_da, &coords);
793: DMDAVecGetArray(cda, coords, &_coords);
795: /* setup for coefficients */
796: DMCreateLocalVector(properties_da, &local_properties);
797: DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties);
798: DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties);
799: DMDAVecGetArray(properties_da, local_properties, &props);
801: DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL);
802: DMDAGetElementsSizes(stokes_da, &mx, &my, NULL);
803: for (ej = sey; ej < sey + my; ej++) {
804: for (ei = sex; ei < sex + mx; ei++) {
805: /* get coords for the element */
806: GetElementCoords(_coords, ei, ej, el_coords);
808: /* get coefficients for the element */
809: prop_eta = props[ej][ei].eta;
811: /* initialise element stiffness matrix */
812: PetscMemzero(Ae, sizeof(Ae));
813: PetscMemzero(Ge, sizeof(Ge));
814: PetscMemzero(De, sizeof(De));
815: PetscMemzero(Ce, sizeof(Ce));
817: /* form element stiffness matrix */
818: FormStressOperatorQ1(Ae, el_coords, prop_eta);
819: FormGradientOperatorQ1(Ge, el_coords);
820: FormScaledMassMatrixOperatorQ1(Ce, el_coords, prop_eta);
822: /* insert element matrix into global matrix */
823: DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej);
824: MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES);
825: MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES);
826: MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES);
827: }
828: }
829: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
830: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
832: DMDAVecRestoreArray(cda, coords, &_coords);
834: DMDAVecRestoreArray(properties_da, local_properties, &props);
835: VecDestroy(&local_properties);
836: return 0;
837: }
839: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F, MatStencil u_eqn[], MatStencil p_eqn[], PetscScalar Fe_u[], PetscScalar Fe_p[])
840: {
841: PetscInt n;
844: for (n = 0; n < 4; n++) {
845: 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];
846: 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];
847: 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];
848: }
849: return 0;
850: }
852: static PetscErrorCode AssembleF_Stokes(Vec F, DM stokes_da, DM properties_da, Vec properties)
853: {
854: DM cda;
855: Vec coords;
856: DMDACoor2d **_coords;
857: MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
858: MatStencil p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
859: PetscInt sex, sey, mx, my;
860: PetscInt ei, ej;
861: PetscScalar Fe[NODES_PER_EL * U_DOFS];
862: PetscScalar He[NODES_PER_EL * P_DOFS];
863: PetscScalar el_coords[NODES_PER_EL * NSD];
864: Vec local_properties;
865: GaussPointCoefficients **props;
866: PetscScalar *prop_fx, *prop_fy;
867: Vec local_F;
868: StokesDOF **ff;
871: /* setup for coords */
872: DMGetCoordinateDM(stokes_da, &cda);
873: DMGetCoordinatesLocal(stokes_da, &coords);
874: DMDAVecGetArray(cda, coords, &_coords);
876: /* setup for coefficients */
877: DMGetLocalVector(properties_da, &local_properties);
878: DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties);
879: DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties);
880: DMDAVecGetArray(properties_da, local_properties, &props);
882: /* get access to the vector */
883: DMGetLocalVector(stokes_da, &local_F);
884: VecZeroEntries(local_F);
885: DMDAVecGetArray(stokes_da, local_F, &ff);
887: DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL);
888: DMDAGetElementsSizes(stokes_da, &mx, &my, NULL);
889: for (ej = sey; ej < sey + my; ej++) {
890: for (ei = sex; ei < sex + mx; ei++) {
891: /* get coords for the element */
892: GetElementCoords(_coords, ei, ej, el_coords);
894: /* get coefficients for the element */
895: prop_fx = props[ej][ei].fx;
896: prop_fy = props[ej][ei].fy;
898: /* initialise element stiffness matrix */
899: PetscMemzero(Fe, sizeof(Fe));
900: PetscMemzero(He, sizeof(He));
902: /* form element stiffness matrix */
903: FormMomentumRhsQ1(Fe, el_coords, prop_fx, prop_fy);
905: /* insert element matrix into global matrix */
906: DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej);
908: DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, p_eqn, Fe, He);
909: }
910: }
912: DMDAVecRestoreArray(stokes_da, local_F, &ff);
913: DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F);
914: DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F);
915: DMRestoreLocalVector(stokes_da, &local_F);
917: DMDAVecRestoreArray(cda, coords, &_coords);
919: DMDAVecRestoreArray(properties_da, local_properties, &props);
920: DMRestoreLocalVector(properties_da, &local_properties);
921: return 0;
922: }
924: static PetscErrorCode DMDACreateSolCx(PetscReal eta0, PetscReal eta1, PetscReal xc, PetscInt nz, PetscInt mx, PetscInt my, DM *_da, Vec *_X)
925: {
926: DM da, cda;
927: Vec X;
928: StokesDOF **_stokes;
929: Vec coords;
930: DMDACoor2d **_coords;
931: PetscInt si, sj, ei, ej, i, j;
934: 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);
935: DMSetFromOptions(da);
936: DMSetUp(da);
937: DMDASetFieldName(da, 0, "anlytic_Vx");
938: DMDASetFieldName(da, 1, "anlytic_Vy");
939: DMDASetFieldName(da, 2, "analytic_P");
941: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0., 0.);
943: DMGetCoordinatesLocal(da, &coords);
944: DMGetCoordinateDM(da, &cda);
945: DMDAVecGetArray(cda, coords, &_coords);
947: DMCreateGlobalVector(da, &X);
948: DMDAVecGetArray(da, X, &_stokes);
950: DMDAGetCorners(da, &si, &sj, 0, &ei, &ej, 0);
951: for (j = sj; j < sj + ej; j++) {
952: for (i = si; i < si + ei; i++) {
953: PetscReal pos[2], pressure, vel[2], total_stress[3], strain_rate[3];
955: pos[0] = PetscRealPart(_coords[j][i].x);
956: pos[1] = PetscRealPart(_coords[j][i].y);
958: evaluate_solCx(pos, eta0, eta1, xc, nz, vel, &pressure, total_stress, strain_rate);
960: _stokes[j][i].u_dof = vel[0];
961: _stokes[j][i].v_dof = vel[1];
962: _stokes[j][i].p_dof = pressure;
963: }
964: }
965: DMDAVecRestoreArray(da, X, &_stokes);
966: DMDAVecRestoreArray(cda, coords, &_coords);
968: *_da = da;
969: *_X = X;
970: return 0;
971: }
973: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields, PetscInt ei, PetscInt ej, StokesDOF nodal_fields[])
974: {
976: /* get the nodal fields */
977: nodal_fields[0].u_dof = fields[ej][ei].u_dof;
978: nodal_fields[0].v_dof = fields[ej][ei].v_dof;
979: nodal_fields[0].p_dof = fields[ej][ei].p_dof;
980: nodal_fields[1].u_dof = fields[ej + 1][ei].u_dof;
981: nodal_fields[1].v_dof = fields[ej + 1][ei].v_dof;
982: nodal_fields[1].p_dof = fields[ej + 1][ei].p_dof;
983: nodal_fields[2].u_dof = fields[ej + 1][ei + 1].u_dof;
984: nodal_fields[2].v_dof = fields[ej + 1][ei + 1].v_dof;
985: nodal_fields[2].p_dof = fields[ej + 1][ei + 1].p_dof;
986: nodal_fields[3].u_dof = fields[ej][ei + 1].u_dof;
987: nodal_fields[3].v_dof = fields[ej][ei + 1].v_dof;
988: nodal_fields[3].p_dof = fields[ej][ei + 1].p_dof;
989: return 0;
990: }
992: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da, Vec X, Vec X_analytic)
993: {
994: DM cda;
995: Vec coords, X_analytic_local, X_local;
996: DMDACoor2d **_coords;
997: PetscInt sex, sey, mx, my;
998: PetscInt ei, ej;
999: PetscScalar el_coords[NODES_PER_EL * NSD];
1000: StokesDOF **stokes_analytic, **stokes;
1001: StokesDOF stokes_analytic_e[4], stokes_e[4];
1003: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
1004: PetscScalar Ni_p[NODES_PER_EL];
1005: PetscInt ngp;
1006: PetscScalar gp_xi[GAUSS_POINTS][2];
1007: PetscScalar gp_weight[GAUSS_POINTS];
1008: PetscInt p, i;
1009: PetscScalar J_p, fac;
1010: PetscScalar h, p_e_L2, u_e_L2, u_e_H1, p_L2, u_L2, u_H1, tp_L2, tu_L2, tu_H1;
1011: PetscInt M;
1012: PetscReal xymin[2], xymax[2];
1015: /* define quadrature rule */
1016: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
1018: /* setup for coords */
1019: DMGetCoordinateDM(stokes_da, &cda);
1020: DMGetCoordinatesLocal(stokes_da, &coords);
1021: DMDAVecGetArray(cda, coords, &_coords);
1023: /* setup for analytic */
1024: DMCreateLocalVector(stokes_da, &X_analytic_local);
1025: DMGlobalToLocalBegin(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local);
1026: DMGlobalToLocalEnd(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local);
1027: DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic);
1029: /* setup for solution */
1030: DMCreateLocalVector(stokes_da, &X_local);
1031: DMGlobalToLocalBegin(stokes_da, X, INSERT_VALUES, X_local);
1032: DMGlobalToLocalEnd(stokes_da, X, INSERT_VALUES, X_local);
1033: DMDAVecGetArray(stokes_da, X_local, &stokes);
1035: DMDAGetInfo(stokes_da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
1036: DMGetBoundingBox(stokes_da, xymin, xymax);
1038: h = (xymax[0] - xymin[0]) / ((PetscReal)M);
1040: tp_L2 = tu_L2 = tu_H1 = 0.0;
1042: DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL);
1043: DMDAGetElementsSizes(stokes_da, &mx, &my, NULL);
1044: for (ej = sey; ej < sey + my; ej++) {
1045: for (ei = sex; ei < sex + mx; ei++) {
1046: /* get coords for the element */
1047: GetElementCoords(_coords, ei, ej, el_coords);
1048: StokesDAGetNodalFields(stokes, ei, ej, stokes_e);
1049: StokesDAGetNodalFields(stokes_analytic, ei, ej, stokes_analytic_e);
1051: /* evaluate integral */
1052: p_e_L2 = 0.0;
1053: u_e_L2 = 0.0;
1054: u_e_H1 = 0.0;
1055: for (p = 0; p < ngp; p++) {
1056: ConstructQ12D_Ni(gp_xi[p], Ni_p);
1057: ConstructQ12D_GNi(gp_xi[p], GNi_p);
1058: ConstructQ12D_GNx(GNi_p, GNx_p, el_coords, &J_p);
1059: fac = gp_weight[p] * J_p;
1061: for (i = 0; i < NODES_PER_EL; i++) {
1062: PetscScalar u_error, v_error;
1064: 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);
1066: u_error = stokes_e[i].u_dof - stokes_analytic_e[i].u_dof;
1067: v_error = stokes_e[i].v_dof - stokes_analytic_e[i].v_dof;
1068: u_e_L2 += fac * Ni_p[i] * (u_error * u_error + v_error * v_error);
1070: u_e_H1 = u_e_H1 + fac * (GNx_p[0][i] * u_error * GNx_p[0][i] * u_error /* du/dx */
1071: + GNx_p[1][i] * u_error * GNx_p[1][i] * u_error /* du/dy */
1072: + GNx_p[0][i] * v_error * GNx_p[0][i] * v_error /* dv/dx */
1073: + GNx_p[1][i] * v_error * GNx_p[1][i] * v_error); /* dv/dy */
1074: }
1075: }
1077: tp_L2 += p_e_L2;
1078: tu_L2 += u_e_L2;
1079: tu_H1 += u_e_H1;
1080: }
1081: }
1082: MPI_Allreduce(&tp_L2, &p_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD);
1083: MPI_Allreduce(&tu_L2, &u_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD);
1084: MPI_Allreduce(&tu_H1, &u_H1, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD);
1085: p_L2 = PetscSqrtScalar(p_L2);
1086: u_L2 = PetscSqrtScalar(u_L2);
1087: u_H1 = PetscSqrtScalar(u_H1);
1089: 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));
1091: DMDAVecRestoreArray(cda, coords, &_coords);
1093: DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic);
1094: VecDestroy(&X_analytic_local);
1095: DMDAVecRestoreArray(stokes_da, X_local, &stokes);
1096: VecDestroy(&X_local);
1097: return 0;
1098: }
1100: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx, PetscInt my)
1101: {
1102: DM da_Stokes, da_prop;
1103: PetscInt u_dof, p_dof, dof, stencil_width;
1104: Mat A, B;
1105: DM prop_cda, vel_cda;
1106: Vec prop_coords, vel_coords;
1107: PetscInt si, sj, nx, ny, i, j, p;
1108: Vec f, X;
1109: PetscInt prop_dof, prop_stencil_width;
1110: Vec properties, l_properties;
1111: PetscReal dx, dy;
1112: PetscInt M, N;
1113: DMDACoor2d **_prop_coords, **_vel_coords;
1114: GaussPointCoefficients **element_props;
1115: PetscInt its;
1116: KSP ksp_S;
1117: PetscInt coefficient_structure = 0;
1118: PetscInt cpu_x, cpu_y, *lx = NULL, *ly = NULL;
1119: PetscBool use_gp_coords = PETSC_FALSE, set, output_gnuplot = PETSC_FALSE, glvis = PETSC_FALSE, change = PETSC_FALSE;
1120: char filename[PETSC_MAX_PATH_LEN];
1124: PetscOptionsGetBool(NULL, NULL, "-gnuplot", &output_gnuplot, NULL);
1125: PetscOptionsGetBool(NULL, NULL, "-glvis", &glvis, NULL);
1126: PetscOptionsGetBool(NULL, NULL, "-change", &change, NULL);
1128: /* Generate the da for velocity and pressure */
1129: /*
1130: We use Q1 elements for the temperature.
1131: FEM has a 9-point stencil (BOX) or connectivity pattern
1132: Num nodes in each direction is mx+1, my+1
1133: */
1134: u_dof = U_DOFS; /* Vx, Vy - velocities */
1135: p_dof = P_DOFS; /* p - pressure */
1136: dof = u_dof + p_dof;
1137: stencil_width = 1;
1138: 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);
1140: DMSetMatType(da_Stokes, MATAIJ);
1141: DMSetFromOptions(da_Stokes);
1142: DMSetUp(da_Stokes);
1143: DMDASetFieldName(da_Stokes, 0, "Vx");
1144: DMDASetFieldName(da_Stokes, 1, "Vy");
1145: DMDASetFieldName(da_Stokes, 2, "P");
1147: /* unit box [0,1] x [0,1] */
1148: DMDASetUniformCoordinates(da_Stokes, 0.0, 1.0, 0.0, 1.0, 0., 0.);
1150: /* Generate element properties, we will assume all material properties are constant over the element */
1151: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1152: DMDAGetInfo(da_Stokes, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0);
1153: DMDAGetElementOwnershipRanges2d(da_Stokes, &lx, &ly);
1155: prop_dof = (int)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
1156: prop_stencil_width = 0;
1157: 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);
1158: DMSetFromOptions(da_prop);
1159: DMSetUp(da_prop);
1160: PetscFree(lx);
1161: PetscFree(ly);
1163: /* define centroid positions */
1164: DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
1165: dx = 1.0 / ((PetscReal)(M));
1166: dy = 1.0 / ((PetscReal)(N));
1168: 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);
1170: /* define coefficients */
1171: PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL);
1173: DMCreateGlobalVector(da_prop, &properties);
1174: DMCreateLocalVector(da_prop, &l_properties);
1175: DMDAVecGetArray(da_prop, l_properties, &element_props);
1177: DMGetCoordinateDM(da_prop, &prop_cda);
1178: DMGetCoordinatesLocal(da_prop, &prop_coords);
1179: DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords);
1181: DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0);
1183: DMGetCoordinateDM(da_Stokes, &vel_cda);
1184: DMGetCoordinatesLocal(da_Stokes, &vel_coords);
1185: DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords);
1187: /* interpolate the coordinates */
1188: for (j = sj; j < sj + ny; j++) {
1189: for (i = si; i < si + nx; i++) {
1190: PetscInt ngp;
1191: PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
1192: PetscScalar el_coords[8];
1194: GetElementCoords(_vel_coords, i, j, el_coords);
1195: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
1197: for (p = 0; p < GAUSS_POINTS; p++) {
1198: PetscScalar gp_x, gp_y;
1199: PetscInt n;
1200: PetscScalar xi_p[2], Ni_p[4];
1202: xi_p[0] = gp_xi[p][0];
1203: xi_p[1] = gp_xi[p][1];
1204: ConstructQ12D_Ni(xi_p, Ni_p);
1206: gp_x = 0.0;
1207: gp_y = 0.0;
1208: for (n = 0; n < NODES_PER_EL; n++) {
1209: gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
1210: gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
1211: }
1212: element_props[j][i].gp_coords[2 * p] = gp_x;
1213: element_props[j][i].gp_coords[2 * p + 1] = gp_y;
1214: }
1215: }
1216: }
1218: /* define the coefficients */
1219: PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, NULL);
1221: for (j = sj; j < sj + ny; j++) {
1222: for (i = si; i < si + nx; i++) {
1223: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1224: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1225: PetscReal coord_x, coord_y;
1227: if (coefficient_structure == 0) {
1228: PetscReal opts_eta0, opts_eta1, opts_xc;
1229: PetscInt opts_nz;
1231: opts_eta0 = 1.0;
1232: opts_eta1 = 1.0;
1233: opts_xc = 0.5;
1234: opts_nz = 1;
1236: PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL);
1237: PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL);
1238: PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL);
1239: PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL);
1241: for (p = 0; p < GAUSS_POINTS; p++) {
1242: coord_x = centroid_x;
1243: coord_y = centroid_y;
1244: if (use_gp_coords) {
1245: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1246: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1247: }
1249: element_props[j][i].eta[p] = opts_eta0;
1250: if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;
1252: element_props[j][i].fx[p] = 0.0;
1253: element_props[j][i].fy[p] = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1254: }
1255: } else if (coefficient_structure == 1) { /* square sinker */
1256: PetscReal opts_eta0, opts_eta1, opts_dx, opts_dy;
1258: opts_eta0 = 1.0;
1259: opts_eta1 = 1.0;
1260: opts_dx = 0.50;
1261: opts_dy = 0.50;
1263: PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL);
1264: PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL);
1265: PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL);
1266: PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, 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: if ((coord_x > -0.5 * opts_dx + 0.5) && (coord_x < 0.5 * opts_dx + 0.5)) {
1281: if ((coord_y > -0.5 * opts_dy + 0.5) && (coord_y < 0.5 * opts_dy + 0.5)) {
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: }
1288: } else if (coefficient_structure == 2) { /* circular sinker */
1289: PetscReal opts_eta0, opts_eta1, opts_r, radius2;
1291: opts_eta0 = 1.0;
1292: opts_eta1 = 1.0;
1293: opts_r = 0.25;
1295: PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL);
1296: PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL);
1297: PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL);
1299: for (p = 0; p < GAUSS_POINTS; p++) {
1300: coord_x = centroid_x;
1301: coord_y = centroid_y;
1302: if (use_gp_coords) {
1303: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1304: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1305: }
1307: element_props[j][i].eta[p] = opts_eta0;
1308: element_props[j][i].fx[p] = 0.0;
1309: element_props[j][i].fy[p] = 0.0;
1311: radius2 = (coord_x - 0.5) * (coord_x - 0.5) + (coord_y - 0.5) * (coord_y - 0.5);
1312: if (radius2 < opts_r * opts_r) {
1313: element_props[j][i].eta[p] = opts_eta1;
1314: element_props[j][i].fx[p] = 0.0;
1315: element_props[j][i].fy[p] = -1.0;
1316: }
1317: }
1318: } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1319: PetscReal opts_eta0, opts_eta1, opts_r, opts_dx, opts_dy, opts_c0x, opts_c0y, opts_s0x, opts_s0y, opts_phi, radius2;
1321: opts_eta0 = 1.0;
1322: opts_eta1 = 1.0;
1323: opts_r = 0.25;
1324: opts_c0x = 0.35; /* circle center */
1325: opts_c0y = 0.35;
1326: opts_s0x = 0.7; /* square center */
1327: opts_s0y = 0.7;
1328: opts_dx = 0.25;
1329: opts_dy = 0.25;
1330: opts_phi = 25;
1332: PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL);
1333: PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL);
1334: PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL);
1335: PetscOptionsGetReal(NULL, NULL, "-sinker_c0x", &opts_c0x, NULL);
1336: PetscOptionsGetReal(NULL, NULL, "-sinker_c0y", &opts_c0y, NULL);
1337: PetscOptionsGetReal(NULL, NULL, "-sinker_s0x", &opts_s0x, NULL);
1338: PetscOptionsGetReal(NULL, NULL, "-sinker_s0y", &opts_s0y, NULL);
1339: PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL);
1340: PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL);
1341: PetscOptionsGetReal(NULL, NULL, "-sinker_phi", &opts_phi, NULL);
1342: opts_phi *= PETSC_PI / 180;
1344: for (p = 0; p < GAUSS_POINTS; p++) {
1345: coord_x = centroid_x;
1346: coord_y = centroid_y;
1347: if (use_gp_coords) {
1348: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1349: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1350: }
1352: element_props[j][i].eta[p] = opts_eta0;
1353: element_props[j][i].fx[p] = 0.0;
1354: element_props[j][i].fy[p] = -0.2;
1356: radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1357: if (radius2 < opts_r * opts_r || (PetscAbs(+(coord_x - opts_s0x) * PetscCosReal(opts_phi) + (coord_y - opts_s0y) * PetscSinReal(opts_phi)) < opts_dx / 2 && PetscAbs(-(coord_x - opts_s0x) * PetscSinReal(opts_phi) + (coord_y - opts_s0y) * PetscCosReal(opts_phi)) < opts_dy / 2)) {
1358: element_props[j][i].eta[p] = opts_eta1;
1359: element_props[j][i].fx[p] = 0.0;
1360: element_props[j][i].fy[p] = -1.0;
1361: }
1362: }
1363: } else if (coefficient_structure == 4) { /* subdomain jump */
1364: PetscReal opts_mag, opts_eta0;
1365: PetscInt opts_nz, px, py;
1366: PetscBool jump;
1368: opts_mag = 1.0;
1369: opts_eta0 = 1.0;
1370: opts_nz = 1;
1372: PetscOptionsGetReal(NULL, NULL, "-jump_eta0", &opts_eta0, NULL);
1373: PetscOptionsGetReal(NULL, NULL, "-jump_magnitude", &opts_mag, NULL);
1374: PetscOptionsGetInt(NULL, NULL, "-jump_nz", &opts_nz, NULL);
1375: DMDAGetInfo(da_Stokes, NULL, NULL, NULL, NULL, &px, &py, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
1376: if (px % 2) {
1377: jump = (PetscBool)(PetscGlobalRank % 2);
1378: } else {
1379: jump = (PetscBool)((PetscGlobalRank / px) % 2 ? PetscGlobalRank % 2 : !(PetscGlobalRank % 2));
1380: }
1381: for (p = 0; p < GAUSS_POINTS; p++) {
1382: coord_x = centroid_x;
1383: coord_y = centroid_y;
1384: if (use_gp_coords) {
1385: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1386: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1387: }
1389: element_props[j][i].eta[p] = jump ? PetscPowReal(10.0, opts_mag) : opts_eta0;
1390: element_props[j][i].fx[p] = 0.0;
1391: element_props[j][i].fy[p] = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1392: }
1393: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
1394: }
1395: }
1396: DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords);
1398: DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords);
1400: DMDAVecRestoreArray(da_prop, l_properties, &element_props);
1401: DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties);
1402: DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties);
1404: if (output_gnuplot) {
1405: DMDACoordViewGnuplot2d(da_Stokes, "mesh");
1406: DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for Stokes eqn.", "properties");
1407: }
1409: if (glvis) {
1410: Vec glv_prop, etaf;
1411: PetscViewer view;
1412: PetscInt dim;
1413: const char *fec = {"FiniteElementCollection: L2_2D_P1"};
1415: DMGetDimension(da_Stokes, &dim);
1416: VecCreateSeq(PETSC_COMM_SELF, GAUSS_POINTS * mx * mx, &etaf);
1417: PetscObjectSetName((PetscObject)etaf, "viscosity");
1418: PetscViewerGLVisOpen(PETSC_COMM_WORLD, PETSC_VIEWER_GLVIS_SOCKET, NULL, PETSC_DECIDE, &view);
1419: PetscViewerGLVisSetFields(view, 1, &fec, &dim, glvis_extract_eta, (PetscObject *)&etaf, da_prop, NULL);
1420: DMGetLocalVector(da_prop, &glv_prop);
1421: DMGlobalToLocalBegin(da_prop, properties, INSERT_VALUES, glv_prop);
1422: DMGlobalToLocalEnd(da_prop, properties, INSERT_VALUES, glv_prop);
1423: VecSetDM(etaf, da_Stokes);
1424: VecView(glv_prop, view);
1425: PetscViewerDestroy(&view);
1426: DMRestoreLocalVector(da_prop, &glv_prop);
1427: VecDestroy(&etaf);
1428: }
1430: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1431: DMCreateMatrix(da_Stokes, &A);
1432: DMCreateMatrix(da_Stokes, &B);
1433: DMCreateGlobalVector(da_Stokes, &f);
1434: DMCreateGlobalVector(da_Stokes, &X);
1436: /* assemble A11 */
1437: MatZeroEntries(A);
1438: MatZeroEntries(B);
1439: VecZeroEntries(f);
1441: AssembleA_Stokes(A, da_Stokes, da_prop, properties);
1442: AssembleA_PCStokes(B, da_Stokes, da_prop, properties);
1443: /* build force vector */
1444: AssembleF_Stokes(f, da_Stokes, da_prop, properties);
1446: DMDABCApplyFreeSlip(da_Stokes, A, f);
1447: DMDABCApplyFreeSlip(da_Stokes, B, NULL);
1449: /* SOLVE */
1450: KSPCreate(PETSC_COMM_WORLD, &ksp_S);
1451: KSPSetOptionsPrefix(ksp_S, "stokes_");
1452: KSPSetDM(ksp_S, da_Stokes);
1453: KSPSetDMActive(ksp_S, PETSC_FALSE);
1454: KSPSetOperators(ksp_S, A, B);
1455: KSPSetFromOptions(ksp_S);
1456: {
1457: PC pc;
1458: const PetscInt ufields[] = {0, 1}, pfields[1] = {2};
1460: KSPGetPC(ksp_S, &pc);
1461: PCFieldSplitSetBlockSize(pc, 3);
1462: PCFieldSplitSetFields(pc, "u", 2, ufields, ufields);
1463: PCFieldSplitSetFields(pc, "p", 1, pfields, pfields);
1464: }
1466: {
1467: PC pc;
1468: PetscBool same = PETSC_FALSE;
1469: KSPGetPC(ksp_S, &pc);
1470: PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &same);
1471: if (same) {
1472: PetscBool usedivmat = PETSC_FALSE;
1473: KSPSetOperators(ksp_S, A, A);
1475: PetscOptionsGetBool(NULL, NULL, "-stokes_pc_bddc_use_divergence", &usedivmat, NULL);
1476: if (usedivmat) {
1477: IS *fields, vel;
1478: PetscInt i, nf;
1480: DMCreateFieldDecomposition(da_Stokes, &nf, NULL, &fields, NULL);
1481: ISConcatenate(PETSC_COMM_WORLD, 2, fields, &vel);
1482: MatZeroRowsIS(B, fields[2], 1.0, NULL, NULL); /* we put 1.0 on the diagonal to pick the pressure average too */
1483: MatTranspose(B, MAT_INPLACE_MATRIX, &B);
1484: MatZeroRowsIS(B, vel, 0.0, NULL, NULL);
1485: ISDestroy(&vel);
1486: PCBDDCSetDivergenceMat(pc, B, PETSC_FALSE, NULL);
1487: for (i = 0; i < nf; i++) ISDestroy(&fields[i]);
1488: PetscFree(fields);
1489: }
1490: }
1491: }
1493: {
1494: PC pc_S;
1495: KSP *sub_ksp, ksp_U;
1496: PetscInt nsplits;
1497: DM da_U;
1498: PetscBool is_pcfs;
1500: KSPSetUp(ksp_S);
1501: KSPGetPC(ksp_S, &pc_S);
1503: is_pcfs = PETSC_FALSE;
1504: PetscObjectTypeCompare((PetscObject)pc_S, PCFIELDSPLIT, &is_pcfs);
1506: if (is_pcfs) {
1507: PCFieldSplitGetSubKSP(pc_S, &nsplits, &sub_ksp);
1508: ksp_U = sub_ksp[0];
1509: PetscFree(sub_ksp);
1511: if (nsplits == 2) {
1512: DMDACreateCompatibleDMDA(da_Stokes, 2, &da_U);
1514: KSPSetDM(ksp_U, da_U);
1515: KSPSetDMActive(ksp_U, PETSC_FALSE);
1516: DMDestroy(&da_U);
1517: if (change) {
1518: const char opt[] = "-stokes_fieldsplit_u_pc_factor_mat_solver_type mumps";
1519: PetscOptionsInsertString(NULL, opt);
1520: KSPSetFromOptions(ksp_U);
1521: }
1522: KSPSetUp(ksp_U);
1523: }
1524: }
1525: }
1527: KSPSolve(ksp_S, f, X);
1529: PetscOptionsGetString(NULL, NULL, "-o", filename, sizeof(filename), &set);
1530: if (set) {
1531: char *ext;
1532: PetscViewer viewer;
1533: PetscBool flg;
1534: PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
1535: PetscStrrchr(filename, '.', &ext);
1536: PetscStrcmp("vts", ext, &flg);
1537: if (flg) {
1538: PetscViewerSetType(viewer, PETSCVIEWERVTK);
1539: } else {
1540: PetscViewerSetType(viewer, PETSCVIEWERBINARY);
1541: }
1542: PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
1543: PetscViewerFileSetName(viewer, filename);
1544: VecView(X, viewer);
1545: PetscViewerDestroy(&viewer);
1546: }
1547: if (output_gnuplot) DMDAViewGnuplot2d(da_Stokes, X, "Velocity solution for Stokes eqn.", "X");
1549: if (glvis) {
1550: PetscViewer view;
1552: PetscViewerCreate(PETSC_COMM_WORLD, &view);
1553: PetscViewerSetType(view, PETSCVIEWERGLVIS);
1554: VecView(X, view);
1555: PetscViewerDestroy(&view);
1556: }
1558: KSPGetIterationNumber(ksp_S, &its);
1560: if (coefficient_structure == 0) {
1561: PetscReal opts_eta0, opts_eta1, opts_xc;
1562: PetscInt opts_nz, N;
1563: DM da_Stokes_analytic;
1564: Vec X_analytic;
1565: PetscReal nrm1[3], nrm2[3], nrmI[3];
1567: opts_eta0 = 1.0;
1568: opts_eta1 = 1.0;
1569: opts_xc = 0.5;
1570: opts_nz = 1;
1572: PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL);
1573: PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL);
1574: PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL);
1575: PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL);
1577: DMDACreateSolCx(opts_eta0, opts_eta1, opts_xc, opts_nz, mx, my, &da_Stokes_analytic, &X_analytic);
1578: if (output_gnuplot) DMDAViewGnuplot2d(da_Stokes_analytic, X_analytic, "Analytic solution for Stokes eqn.", "X_analytic");
1579: DMDAIntegrateErrors(da_Stokes_analytic, X, X_analytic);
1581: VecAXPY(X_analytic, -1.0, X);
1582: VecGetSize(X_analytic, &N);
1583: N = N / 3;
1585: VecStrideNorm(X_analytic, 0, NORM_1, &nrm1[0]);
1586: VecStrideNorm(X_analytic, 0, NORM_2, &nrm2[0]);
1587: VecStrideNorm(X_analytic, 0, NORM_INFINITY, &nrmI[0]);
1589: VecStrideNorm(X_analytic, 1, NORM_1, &nrm1[1]);
1590: VecStrideNorm(X_analytic, 1, NORM_2, &nrm2[1]);
1591: VecStrideNorm(X_analytic, 1, NORM_INFINITY, &nrmI[1]);
1593: VecStrideNorm(X_analytic, 2, NORM_1, &nrm1[2]);
1594: VecStrideNorm(X_analytic, 2, NORM_2, &nrm2[2]);
1595: VecStrideNorm(X_analytic, 2, NORM_INFINITY, &nrmI[2]);
1597: DMDestroy(&da_Stokes_analytic);
1598: VecDestroy(&X_analytic);
1599: }
1601: KSPDestroy(&ksp_S);
1602: VecDestroy(&X);
1603: VecDestroy(&f);
1604: MatDestroy(&A);
1605: MatDestroy(&B);
1607: DMDestroy(&da_Stokes);
1608: DMDestroy(&da_prop);
1610: VecDestroy(&properties);
1611: VecDestroy(&l_properties);
1612: return 0;
1613: }
1615: int main(int argc, char **args)
1616: {
1617: PetscInt mx, my;
1620: PetscInitialize(&argc, &args, (char *)0, help);
1621: mx = my = 10;
1622: PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL);
1623: PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL);
1624: solve_stokes_2d_coupled(mx, my);
1625: PetscFinalize();
1626: return 0;
1627: }
1629: /* -------------------------- helpers for boundary conditions -------------------------------- */
1630: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
1631: {
1632: DM cda;
1633: Vec coords;
1634: PetscInt si, sj, nx, ny, i, j;
1635: PetscInt M, N;
1636: DMDACoor2d **_coords;
1637: const PetscInt *g_idx;
1638: PetscInt *bc_global_ids;
1639: PetscScalar *bc_vals;
1640: PetscInt nbcs;
1641: PetscInt n_dofs;
1642: ISLocalToGlobalMapping ltogm;
1645: DMGetLocalToGlobalMapping(da, <ogm);
1646: ISLocalToGlobalMappingGetIndices(ltogm, &g_idx);
1648: DMGetCoordinateDM(da, &cda);
1649: DMGetCoordinatesLocal(da, &coords);
1650: DMDAVecGetArray(cda, coords, &_coords);
1651: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
1652: DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);
1654: PetscMalloc1(ny * n_dofs, &bc_global_ids);
1655: PetscMalloc1(ny * n_dofs, &bc_vals);
1657: /* init the entries to -1 so VecSetValues will ignore them */
1658: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1660: i = nx - 1;
1661: for (j = 0; j < ny; j++) {
1662: PetscInt local_id;
1664: local_id = i + j * nx;
1666: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1668: bc_vals[j] = 0.0;
1669: }
1670: ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1671: nbcs = 0;
1672: if ((si + nx) == (M)) nbcs = ny;
1674: if (b) {
1675: VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES);
1676: VecAssemblyBegin(b);
1677: VecAssemblyEnd(b);
1678: }
1679: if (A) MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0);
1681: PetscFree(bc_vals);
1682: PetscFree(bc_global_ids);
1684: DMDAVecRestoreArray(cda, coords, &_coords);
1685: return 0;
1686: }
1688: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
1689: {
1690: DM cda;
1691: Vec coords;
1692: PetscInt si, sj, nx, ny, i, j;
1693: PetscInt M, N;
1694: DMDACoor2d **_coords;
1695: const PetscInt *g_idx;
1696: PetscInt *bc_global_ids;
1697: PetscScalar *bc_vals;
1698: PetscInt nbcs;
1699: PetscInt n_dofs;
1700: ISLocalToGlobalMapping ltogm;
1703: DMGetLocalToGlobalMapping(da, <ogm);
1704: ISLocalToGlobalMappingGetIndices(ltogm, &g_idx);
1706: DMGetCoordinateDM(da, &cda);
1707: DMGetCoordinatesLocal(da, &coords);
1708: DMDAVecGetArray(cda, coords, &_coords);
1709: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
1710: DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);
1712: PetscMalloc1(ny * n_dofs, &bc_global_ids);
1713: PetscMalloc1(ny * n_dofs, &bc_vals);
1715: /* init the entries to -1 so VecSetValues will ignore them */
1716: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1718: i = 0;
1719: for (j = 0; j < ny; j++) {
1720: PetscInt local_id;
1722: local_id = i + j * nx;
1724: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1726: bc_vals[j] = 0.0;
1727: }
1728: ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1729: nbcs = 0;
1730: if (si == 0) nbcs = ny;
1732: if (b) {
1733: VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES);
1734: VecAssemblyBegin(b);
1735: VecAssemblyEnd(b);
1736: }
1738: if (A) MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0);
1740: PetscFree(bc_vals);
1741: PetscFree(bc_global_ids);
1743: DMDAVecRestoreArray(cda, coords, &_coords);
1744: return 0;
1745: }
1747: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
1748: {
1749: DM cda;
1750: Vec coords;
1751: PetscInt si, sj, nx, ny, i, j;
1752: PetscInt M, N;
1753: DMDACoor2d **_coords;
1754: const PetscInt *g_idx;
1755: PetscInt *bc_global_ids;
1756: PetscScalar *bc_vals;
1757: PetscInt nbcs;
1758: PetscInt n_dofs;
1759: ISLocalToGlobalMapping ltogm;
1762: DMGetLocalToGlobalMapping(da, <ogm);
1763: ISLocalToGlobalMappingGetIndices(ltogm, &g_idx);
1765: DMGetCoordinateDM(da, &cda);
1766: DMGetCoordinatesLocal(da, &coords);
1767: DMDAVecGetArray(cda, coords, &_coords);
1768: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
1769: DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);
1771: PetscMalloc1(nx, &bc_global_ids);
1772: PetscMalloc1(nx, &bc_vals);
1774: /* init the entries to -1 so VecSetValues will ignore them */
1775: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1777: j = ny - 1;
1778: for (i = 0; i < nx; i++) {
1779: PetscInt local_id;
1781: local_id = i + j * nx;
1783: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1785: bc_vals[i] = 0.0;
1786: }
1787: ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1788: nbcs = 0;
1789: if ((sj + ny) == (N)) nbcs = nx;
1791: if (b) {
1792: VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES);
1793: VecAssemblyBegin(b);
1794: VecAssemblyEnd(b);
1795: }
1796: if (A) MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0);
1798: PetscFree(bc_vals);
1799: PetscFree(bc_global_ids);
1801: DMDAVecRestoreArray(cda, coords, &_coords);
1802: return 0;
1803: }
1805: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
1806: {
1807: DM cda;
1808: Vec coords;
1809: PetscInt si, sj, nx, ny, i, j;
1810: PetscInt M, N;
1811: DMDACoor2d **_coords;
1812: const PetscInt *g_idx;
1813: PetscInt *bc_global_ids;
1814: PetscScalar *bc_vals;
1815: PetscInt nbcs;
1816: PetscInt n_dofs;
1817: ISLocalToGlobalMapping ltogm;
1820: DMGetLocalToGlobalMapping(da, <ogm);
1821: ISLocalToGlobalMappingGetIndices(ltogm, &g_idx);
1823: DMGetCoordinateDM(da, &cda);
1824: DMGetCoordinatesLocal(da, &coords);
1825: DMDAVecGetArray(cda, coords, &_coords);
1826: DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0);
1827: DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0);
1829: PetscMalloc1(nx, &bc_global_ids);
1830: PetscMalloc1(nx, &bc_vals);
1832: /* init the entries to -1 so VecSetValues will ignore them */
1833: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1835: j = 0;
1836: for (i = 0; i < nx; i++) {
1837: PetscInt local_id;
1839: local_id = i + j * nx;
1841: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1843: bc_vals[i] = 0.0;
1844: }
1845: ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1846: nbcs = 0;
1847: if (sj == 0) nbcs = nx;
1849: if (b) {
1850: VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES);
1851: VecAssemblyBegin(b);
1852: VecAssemblyEnd(b);
1853: }
1854: if (A) MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0);
1856: PetscFree(bc_vals);
1857: PetscFree(bc_global_ids);
1859: DMDAVecRestoreArray(cda, coords, &_coords);
1860: return 0;
1861: }
1863: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
1864: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes, Mat A, Vec f)
1865: {
1867: BCApplyZero_NORTH(da_Stokes, 1, A, f);
1868: BCApplyZero_EAST(da_Stokes, 0, A, f);
1869: BCApplyZero_SOUTH(da_Stokes, 1, A, f);
1870: BCApplyZero_WEST(da_Stokes, 0, A, f);
1871: return 0;
1872: }
1874: /*TEST
1876: build:
1877: requires: !complex !single
1879: test:
1880: 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
1882: testset:
1883: 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
1884: test:
1885: suffix: 2
1886: args:
1887: output_file: output/ex43_1.out
1888: test:
1889: requires: mumps
1890: suffix: 2_mumps
1891: args: -change true -stokes_ksp_view
1892: output_file: output/ex43_2_mumps.out
1893: # mumps INFO,INFOG,RINFO,RINFOG may vary on different archs, so keep just a stable one
1894: filter: grep -E -v "(INFOG\([^7]|INFO\(|\[0\])"
1896: test:
1897: suffix: 3
1898: nsize: 4
1899: 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
1901: test:
1902: suffix: 4
1903: nsize: 4
1904: 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
1906: test:
1907: suffix: 5
1908: nsize: 4
1909: 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
1911: test:
1912: suffix: 6
1913: nsize: 8
1914: 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
1916: test:
1917: suffix: bjacobi
1918: nsize: 4
1919: args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type aij -stokes_ksp_converged_reason
1921: test:
1922: suffix: bjacobi_baij
1923: nsize: 4
1924: args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type baij -stokes_ksp_converged_reason
1925: output_file: output/ex43_bjacobi.out
1927: test:
1928: suffix: nested_gmg
1929: nsize: 4
1930: 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
1932: test:
1933: suffix: fetidp
1934: nsize: 8
1935: 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
1937: test:
1938: suffix: fetidp_unsym
1939: nsize: 8
1940: 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
1942: test:
1943: suffix: bddc_stokes_deluxe
1944: nsize: 8
1945: 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
1947: test:
1948: suffix: bddc_stokes_subdomainjump_deluxe
1949: nsize: 9
1950: 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
1952: TEST*/