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;
92: PetscFunctionBeginUser;
93: PetscCall(VecGetDM(Vf[0], &stokes_da));
94: PetscCall(DMDAVecGetArrayRead(properties_da, V, &props));
95: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
96: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
97: PetscCall(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: PetscCall(VecRestoreArray(Vf[0], &array));
105: PetscCall(DMDAVecRestoreArrayRead(properties_da, V, &props));
106: PetscFunctionReturn(PETSC_SUCCESS);
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: {
198: PetscFunctionBeginUser;
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: PetscFunctionReturn(PETSC_SUCCESS);
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;
260: PetscFunctionBeginUser;
261: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
263: PetscCall(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: PetscCall(PetscMalloc1(cpu_x, &LX));
269: PetscCall(PetscMalloc1(cpu_y, &LY));
271: PetscCall(DMDAGetElementsSizes(da, &local_mx, &local_my, NULL));
272: PetscCall(VecCreate(PETSC_COMM_WORLD, &vlx));
273: PetscCall(VecSetSizes(vlx, PETSC_DECIDE, cpu_x));
274: PetscCall(VecSetFromOptions(vlx));
276: PetscCall(VecCreate(PETSC_COMM_WORLD, &vly));
277: PetscCall(VecSetSizes(vly, PETSC_DECIDE, cpu_y));
278: PetscCall(VecSetFromOptions(vly));
280: PetscCall(VecSetValue(vlx, proc_I, (PetscScalar)(local_mx + 1.0e-9), INSERT_VALUES));
281: PetscCall(VecSetValue(vly, proc_J, (PetscScalar)(local_my + 1.0e-9), INSERT_VALUES));
282: PetscCall(VecAssemblyBegin(vlx));
283: PetscCall(VecAssemblyEnd(vlx));
284: PetscCall(VecAssemblyBegin(vly));
285: PetscCall(VecAssemblyEnd(vly));
287: PetscCall(VecScatterCreateToAll(vlx, &ctx, &V_SEQ));
288: PetscCall(VecScatterBegin(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
289: PetscCall(VecScatterEnd(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
290: PetscCall(VecGetArray(V_SEQ, &_a));
291: for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
292: PetscCall(VecRestoreArray(V_SEQ, &_a));
293: PetscCall(VecScatterDestroy(&ctx));
294: PetscCall(VecDestroy(&V_SEQ));
296: PetscCall(VecScatterCreateToAll(vly, &ctx, &V_SEQ));
297: PetscCall(VecScatterBegin(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
298: PetscCall(VecScatterEnd(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
299: PetscCall(VecGetArray(V_SEQ, &_a));
300: for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
301: PetscCall(VecRestoreArray(V_SEQ, &_a));
302: PetscCall(VecScatterDestroy(&ctx));
303: PetscCall(VecDestroy(&V_SEQ));
305: *_lx = LX;
306: *_ly = LY;
308: PetscCall(VecDestroy(&vlx));
309: PetscCall(VecDestroy(&vly));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode DMDACoordViewGnuplot2d(DM da, const char prefix[])
314: {
315: DM cda;
316: Vec coords;
317: DMDACoor2d **_coords;
318: PetscInt si, sj, nx, ny, i, j;
319: FILE *fp;
320: char fname[PETSC_MAX_PATH_LEN];
321: PetscMPIInt rank;
323: PetscFunctionBeginUser;
324: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
325: PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
326: PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
327: PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
329: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### Element geometry for processor %1.4d ### \n", rank));
331: PetscCall(DMGetCoordinateDM(da, &cda));
332: PetscCall(DMGetCoordinatesLocal(da, &coords));
333: PetscCall(DMDAVecGetArrayRead(cda, coords, &_coords));
334: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
335: for (j = sj; j < sj + ny - 1; j++) {
336: for (i = si; i < si + nx - 1; i++) {
337: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
338: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i].x), (double)PetscRealPart(_coords[j + 1][i].y)));
339: PetscCall(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)));
340: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i + 1].x), (double)PetscRealPart(_coords[j][i + 1].y)));
341: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n\n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
342: }
343: }
344: PetscCall(DMDAVecRestoreArrayRead(cda, coords, &_coords));
346: PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: static PetscErrorCode DMDAViewGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
351: {
352: DM cda;
353: Vec coords, local_fields;
354: DMDACoor2d **_coords;
355: FILE *fp;
356: char fname[PETSC_MAX_PATH_LEN];
357: PetscMPIInt rank;
358: PetscInt si, sj, nx, ny, i, j;
359: PetscInt n_dofs, d;
360: PetscScalar *_fields;
362: PetscFunctionBeginUser;
363: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
364: PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
365: PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
366: PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
368: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
369: PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
370: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
371: for (d = 0; d < n_dofs; d++) {
372: const char *field_name;
373: PetscCall(DMDAGetFieldName(da, d, &field_name));
374: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
375: }
376: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));
378: PetscCall(DMGetCoordinateDM(da, &cda));
379: PetscCall(DMGetCoordinatesLocal(da, &coords));
380: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
381: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
383: PetscCall(DMCreateLocalVector(da, &local_fields));
384: PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
385: PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
386: PetscCall(VecGetArray(local_fields, &_fields));
388: for (j = sj; j < sj + ny; j++) {
389: for (i = si; i < si + nx; i++) {
390: PetscScalar coord_x, coord_y;
391: PetscScalar field_d;
393: coord_x = _coords[j][i].x;
394: coord_y = _coords[j][i].y;
396: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
397: for (d = 0; d < n_dofs; d++) {
398: field_d = _fields[n_dofs * ((i - si) + (j - sj) * (nx)) + d];
399: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e ", (double)PetscRealPart(field_d)));
400: }
401: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
402: }
403: }
404: PetscCall(VecRestoreArray(local_fields, &_fields));
405: PetscCall(VecDestroy(&local_fields));
407: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
409: PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
414: {
415: DM cda;
416: Vec local_fields;
417: FILE *fp;
418: char fname[PETSC_MAX_PATH_LEN];
419: PetscMPIInt rank;
420: PetscInt si, sj, nx, ny, i, j, p;
421: PetscInt n_dofs, d;
422: GaussPointCoefficients **_coefficients;
424: PetscFunctionBeginUser;
425: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
426: PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
427: PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
428: PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
430: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
431: PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
432: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
433: for (d = 0; d < n_dofs; d++) {
434: const char *field_name;
435: PetscCall(DMDAGetFieldName(da, d, &field_name));
436: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
437: }
438: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));
440: PetscCall(DMGetCoordinateDM(da, &cda));
441: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
443: PetscCall(DMCreateLocalVector(da, &local_fields));
444: PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
445: PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
446: PetscCall(DMDAVecGetArray(da, local_fields, &_coefficients));
448: for (j = sj; j < sj + ny; j++) {
449: for (i = si; i < si + nx; i++) {
450: PetscScalar coord_x, coord_y;
452: for (p = 0; p < GAUSS_POINTS; p++) {
453: coord_x = _coefficients[j][i].gp_coords[2 * p];
454: coord_y = _coefficients[j][i].gp_coords[2 * p + 1];
456: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
458: PetscCall(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])));
459: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
460: }
461: }
462: }
463: PetscCall(DMDAVecRestoreArray(da, local_fields, &_coefficients));
464: PetscCall(VecDestroy(&local_fields));
466: PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: 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)
471: {
472: PetscInt ij;
473: PetscInt r, c, nc;
475: nc = u_NPE * u_dof;
476: r = w_dof * wi + wd;
477: c = u_dof * ui + ud;
478: ij = r * nc + c;
479: return ij;
480: }
482: /*
483: D = [ 2.eta 0 0 ]
484: [ 0 2.eta 0 ]
485: [ 0 0 eta ]
487: B = [ d_dx 0 ]
488: [ 0 d_dy ]
489: [ d_dy d_dx ]
490: */
491: static void FormStressOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
492: {
493: PetscInt ngp;
494: PetscScalar gp_xi[GAUSS_POINTS][2];
495: PetscScalar gp_weight[GAUSS_POINTS];
496: PetscInt p, i, j, k;
497: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
498: PetscScalar J_p, tildeD[3];
499: PetscScalar B[3][U_DOFS * NODES_PER_EL];
501: /* define quadrature rule */
502: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
504: /* evaluate integral */
505: for (p = 0; p < ngp; p++) {
506: ConstructQ12D_GNi(gp_xi[p], GNi_p);
507: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
509: for (i = 0; i < NODES_PER_EL; i++) {
510: PetscScalar d_dx_i = GNx_p[0][i];
511: PetscScalar d_dy_i = GNx_p[1][i];
513: B[0][2 * i] = d_dx_i;
514: B[0][2 * i + 1] = 0.0;
515: B[1][2 * i] = 0.0;
516: B[1][2 * i + 1] = d_dy_i;
517: B[2][2 * i] = d_dy_i;
518: B[2][2 * i + 1] = d_dx_i;
519: }
521: tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
522: tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
523: tildeD[2] = gp_weight[p] * J_p * eta[p];
525: /* form Bt tildeD B */
526: /*
527: Ke_ij = Bt_ik . D_kl . B_lj
528: = B_ki . D_kl . B_lj
529: = B_ki . D_kk . B_kj
530: */
531: for (i = 0; i < 8; i++) {
532: for (j = 0; j < 8; j++) {
533: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
534: Ke[i + 8 * j] = Ke[i + 8 * j] + B[k][i] * tildeD[k] * B[k][j];
535: }
536: }
537: }
538: }
539: }
541: static void FormGradientOperatorQ1(PetscScalar Ke[], PetscScalar coords[])
542: {
543: PetscInt ngp;
544: PetscScalar gp_xi[GAUSS_POINTS][2];
545: PetscScalar gp_weight[GAUSS_POINTS];
546: PetscInt p, i, j, di;
547: PetscScalar Ni_p[NODES_PER_EL];
548: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
549: PetscScalar J_p, fac;
551: /* define quadrature rule */
552: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
554: /* evaluate integral */
555: for (p = 0; p < ngp; p++) {
556: ConstructQ12D_Ni(gp_xi[p], Ni_p);
557: ConstructQ12D_GNi(gp_xi[p], GNi_p);
558: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
559: fac = gp_weight[p] * J_p;
561: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
562: for (di = 0; di < NSD; di++) { /* u dofs */
563: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
564: PetscInt IJ;
565: IJ = ASS_MAP_wIwDI_uJuDJ(i, di, NODES_PER_EL, 2, j, 0, NODES_PER_EL, 1);
567: Ke[IJ] = Ke[IJ] - GNx_p[di][i] * Ni_p[j] * fac;
568: }
569: }
570: }
571: }
572: }
574: static void FormDivergenceOperatorQ1(PetscScalar De[], PetscScalar coords[])
575: {
576: PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
577: PetscInt i, j;
578: PetscInt nr_g, nc_g;
580: PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));
581: FormGradientOperatorQ1(Ge, coords);
583: nr_g = U_DOFS * NODES_PER_EL;
584: nc_g = P_DOFS * NODES_PER_EL;
586: for (i = 0; i < nr_g; i++) {
587: for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
588: }
589: }
591: static void FormStabilisationOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
592: {
593: PetscInt ngp;
594: PetscScalar gp_xi[GAUSS_POINTS][2];
595: PetscScalar gp_weight[GAUSS_POINTS];
596: PetscInt p, i, j;
597: PetscScalar Ni_p[NODES_PER_EL];
598: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
599: PetscScalar J_p, fac, eta_avg;
601: /* define quadrature rule */
602: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
604: /* evaluate integral */
605: for (p = 0; p < ngp; p++) {
606: ConstructQ12D_Ni(gp_xi[p], Ni_p);
607: ConstructQ12D_GNi(gp_xi[p], GNi_p);
608: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
609: fac = gp_weight[p] * J_p;
611: for (i = 0; i < NODES_PER_EL; i++) {
612: 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);
613: }
614: }
616: /* scale */
617: eta_avg = 0.0;
618: for (p = 0; p < ngp; p++) eta_avg += eta[p];
619: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
620: fac = 1.0 / eta_avg;
621: for (i = 0; i < NODES_PER_EL; i++) {
622: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
623: }
624: }
626: static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
627: {
628: PetscInt ngp;
629: PetscScalar gp_xi[GAUSS_POINTS][2];
630: PetscScalar gp_weight[GAUSS_POINTS];
631: PetscInt p, i, j;
632: PetscScalar Ni_p[NODES_PER_EL];
633: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
634: PetscScalar J_p, fac, eta_avg;
636: /* define quadrature rule */
637: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
639: /* evaluate integral */
640: for (p = 0; p < ngp; p++) {
641: ConstructQ12D_Ni(gp_xi[p], Ni_p);
642: ConstructQ12D_GNi(gp_xi[p], GNi_p);
643: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
644: fac = gp_weight[p] * J_p;
646: for (i = 0; i < NODES_PER_EL; i++) {
647: 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];
648: }
649: }
651: /* scale */
652: eta_avg = 0.0;
653: for (p = 0; p < ngp; p++) eta_avg += eta[p];
654: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
655: fac = 1.0 / eta_avg;
656: for (i = 0; i < NODES_PER_EL; i++) {
657: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
658: }
659: }
661: static void FormMomentumRhsQ1(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
662: {
663: PetscInt ngp;
664: PetscScalar gp_xi[GAUSS_POINTS][2];
665: PetscScalar gp_weight[GAUSS_POINTS];
666: PetscInt p, i;
667: PetscScalar Ni_p[NODES_PER_EL];
668: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
669: PetscScalar J_p, fac;
671: /* define quadrature rule */
672: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
674: /* evaluate integral */
675: for (p = 0; p < ngp; p++) {
676: ConstructQ12D_Ni(gp_xi[p], Ni_p);
677: ConstructQ12D_GNi(gp_xi[p], GNi_p);
678: ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
679: fac = gp_weight[p] * J_p;
681: for (i = 0; i < NODES_PER_EL; i++) {
682: Fe[NSD * i] += fac * Ni_p[i] * fx[p];
683: Fe[NSD * i + 1] += fac * Ni_p[i] * fy[p];
684: }
685: }
686: }
688: static PetscErrorCode GetElementCoords(DMDACoor2d **_coords, PetscInt ei, PetscInt ej, PetscScalar el_coords[])
689: {
690: PetscFunctionBeginUser;
691: /* get coords for the element */
692: el_coords[NSD * 0] = _coords[ej][ei].x;
693: el_coords[NSD * 0 + 1] = _coords[ej][ei].y;
694: el_coords[NSD * 1] = _coords[ej + 1][ei].x;
695: el_coords[NSD * 1 + 1] = _coords[ej + 1][ei].y;
696: el_coords[NSD * 2] = _coords[ej + 1][ei + 1].x;
697: el_coords[NSD * 2 + 1] = _coords[ej + 1][ei + 1].y;
698: el_coords[NSD * 3] = _coords[ej][ei + 1].x;
699: el_coords[NSD * 3 + 1] = _coords[ej][ei + 1].y;
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: static PetscErrorCode AssembleA_Stokes(Mat A, DM stokes_da, DM properties_da, Vec properties)
704: {
705: DM cda;
706: Vec coords;
707: DMDACoor2d **_coords;
708: MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
709: MatStencil p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
710: PetscInt sex, sey, mx, my;
711: PetscInt ei, ej;
712: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
713: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
714: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
715: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
716: PetscScalar el_coords[NODES_PER_EL * NSD];
717: Vec local_properties;
718: GaussPointCoefficients **props;
719: PetscScalar *prop_eta;
721: PetscFunctionBeginUser;
722: /* setup for coords */
723: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
724: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
725: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
727: /* setup for coefficients */
728: PetscCall(DMCreateLocalVector(properties_da, &local_properties));
729: PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
730: PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
731: PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));
733: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
734: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
735: for (ej = sey; ej < sey + my; ej++) {
736: for (ei = sex; ei < sex + mx; ei++) {
737: /* get coords for the element */
738: PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
740: /* get coefficients for the element */
741: prop_eta = props[ej][ei].eta;
743: /* initialise element stiffness matrix */
744: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
745: PetscCall(PetscMemzero(Ge, sizeof(Ge)));
746: PetscCall(PetscMemzero(De, sizeof(De)));
747: PetscCall(PetscMemzero(Ce, sizeof(Ce)));
749: /* form element stiffness matrix */
750: FormStressOperatorQ1(Ae, el_coords, prop_eta);
751: FormGradientOperatorQ1(Ge, el_coords);
752: FormDivergenceOperatorQ1(De, el_coords);
753: FormStabilisationOperatorQ1(Ce, el_coords, prop_eta);
755: /* insert element matrix into global matrix */
756: PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
757: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
758: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
759: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
760: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
761: }
762: }
763: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
764: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
766: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
768: PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
769: PetscCall(VecDestroy(&local_properties));
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: static PetscErrorCode AssembleA_PCStokes(Mat A, DM stokes_da, DM properties_da, Vec properties)
774: {
775: DM cda;
776: Vec coords;
777: DMDACoor2d **_coords;
778: MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
779: MatStencil p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
780: PetscInt sex, sey, mx, my;
781: PetscInt ei, ej;
782: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
783: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
784: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
785: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
786: PetscScalar el_coords[NODES_PER_EL * NSD];
787: Vec local_properties;
788: GaussPointCoefficients **props;
789: PetscScalar *prop_eta;
791: PetscFunctionBeginUser;
792: /* setup for coords */
793: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
794: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
795: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
797: /* setup for coefficients */
798: PetscCall(DMCreateLocalVector(properties_da, &local_properties));
799: PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
800: PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
801: PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));
803: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
804: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
805: for (ej = sey; ej < sey + my; ej++) {
806: for (ei = sex; ei < sex + mx; ei++) {
807: /* get coords for the element */
808: PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
810: /* get coefficients for the element */
811: prop_eta = props[ej][ei].eta;
813: /* initialise element stiffness matrix */
814: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
815: PetscCall(PetscMemzero(Ge, sizeof(Ge)));
816: PetscCall(PetscMemzero(De, sizeof(De)));
817: PetscCall(PetscMemzero(Ce, sizeof(Ce)));
819: /* form element stiffness matrix */
820: FormStressOperatorQ1(Ae, el_coords, prop_eta);
821: FormGradientOperatorQ1(Ge, el_coords);
822: FormScaledMassMatrixOperatorQ1(Ce, el_coords, prop_eta);
824: /* insert element matrix into global matrix */
825: PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
826: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
827: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
828: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
829: }
830: }
831: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
832: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
834: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
836: PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
837: PetscCall(VecDestroy(&local_properties));
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F, MatStencil u_eqn[], MatStencil p_eqn[], PetscScalar Fe_u[], PetscScalar Fe_p[])
842: {
843: PetscInt n;
845: PetscFunctionBeginUser;
846: for (n = 0; n < 4; n++) {
847: 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];
848: 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];
849: 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];
850: }
851: PetscFunctionReturn(PETSC_SUCCESS);
852: }
854: static PetscErrorCode AssembleF_Stokes(Vec F, DM stokes_da, DM properties_da, Vec properties)
855: {
856: DM cda;
857: Vec coords;
858: DMDACoor2d **_coords;
859: MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
860: MatStencil p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
861: PetscInt sex, sey, mx, my;
862: PetscInt ei, ej;
863: PetscScalar Fe[NODES_PER_EL * U_DOFS];
864: PetscScalar He[NODES_PER_EL * P_DOFS];
865: PetscScalar el_coords[NODES_PER_EL * NSD];
866: Vec local_properties;
867: GaussPointCoefficients **props;
868: PetscScalar *prop_fx, *prop_fy;
869: Vec local_F;
870: StokesDOF **ff;
872: PetscFunctionBeginUser;
873: /* setup for coords */
874: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
875: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
876: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
878: /* setup for coefficients */
879: PetscCall(DMGetLocalVector(properties_da, &local_properties));
880: PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
881: PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
882: PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));
884: /* get access to the vector */
885: PetscCall(DMGetLocalVector(stokes_da, &local_F));
886: PetscCall(VecZeroEntries(local_F));
887: PetscCall(DMDAVecGetArray(stokes_da, local_F, &ff));
889: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
890: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
891: for (ej = sey; ej < sey + my; ej++) {
892: for (ei = sex; ei < sex + mx; ei++) {
893: /* get coords for the element */
894: PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
896: /* get coefficients for the element */
897: prop_fx = props[ej][ei].fx;
898: prop_fy = props[ej][ei].fy;
900: /* initialise element stiffness matrix */
901: PetscCall(PetscMemzero(Fe, sizeof(Fe)));
902: PetscCall(PetscMemzero(He, sizeof(He)));
904: /* form element stiffness matrix */
905: FormMomentumRhsQ1(Fe, el_coords, prop_fx, prop_fy);
907: /* insert element matrix into global matrix */
908: PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
910: PetscCall(DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, p_eqn, Fe, He));
911: }
912: }
914: PetscCall(DMDAVecRestoreArray(stokes_da, local_F, &ff));
915: PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
916: PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
917: PetscCall(DMRestoreLocalVector(stokes_da, &local_F));
919: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
921: PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
922: PetscCall(DMRestoreLocalVector(properties_da, &local_properties));
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: static PetscErrorCode DMDACreateSolCx(PetscReal eta0, PetscReal eta1, PetscReal xc, PetscInt nz, PetscInt mx, PetscInt my, DM *_da, Vec *_X)
927: {
928: DM da, cda;
929: Vec X;
930: StokesDOF **_stokes;
931: Vec coords;
932: DMDACoor2d **_coords;
933: PetscInt si, sj, ei, ej, i, j;
935: PetscFunctionBeginUser;
936: PetscCall(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));
937: PetscCall(DMSetFromOptions(da));
938: PetscCall(DMSetUp(da));
939: PetscCall(DMDASetFieldName(da, 0, "anlytic_Vx"));
940: PetscCall(DMDASetFieldName(da, 1, "anlytic_Vy"));
941: PetscCall(DMDASetFieldName(da, 2, "analytic_P"));
943: PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0., 0.));
945: PetscCall(DMGetCoordinatesLocal(da, &coords));
946: PetscCall(DMGetCoordinateDM(da, &cda));
947: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
949: PetscCall(DMCreateGlobalVector(da, &X));
950: PetscCall(DMDAVecGetArray(da, X, &_stokes));
952: PetscCall(DMDAGetCorners(da, &si, &sj, 0, &ei, &ej, 0));
953: for (j = sj; j < sj + ej; j++) {
954: for (i = si; i < si + ei; i++) {
955: PetscReal pos[2], pressure, vel[2], total_stress[3], strain_rate[3];
957: pos[0] = PetscRealPart(_coords[j][i].x);
958: pos[1] = PetscRealPart(_coords[j][i].y);
960: evaluate_solCx(pos, eta0, eta1, xc, nz, vel, &pressure, total_stress, strain_rate);
962: _stokes[j][i].u_dof = vel[0];
963: _stokes[j][i].v_dof = vel[1];
964: _stokes[j][i].p_dof = pressure;
965: }
966: }
967: PetscCall(DMDAVecRestoreArray(da, X, &_stokes));
968: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
970: *_da = da;
971: *_X = X;
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields, PetscInt ei, PetscInt ej, StokesDOF nodal_fields[])
976: {
977: PetscFunctionBeginUser;
978: /* get the nodal fields */
979: nodal_fields[0].u_dof = fields[ej][ei].u_dof;
980: nodal_fields[0].v_dof = fields[ej][ei].v_dof;
981: nodal_fields[0].p_dof = fields[ej][ei].p_dof;
982: nodal_fields[1].u_dof = fields[ej + 1][ei].u_dof;
983: nodal_fields[1].v_dof = fields[ej + 1][ei].v_dof;
984: nodal_fields[1].p_dof = fields[ej + 1][ei].p_dof;
985: nodal_fields[2].u_dof = fields[ej + 1][ei + 1].u_dof;
986: nodal_fields[2].v_dof = fields[ej + 1][ei + 1].v_dof;
987: nodal_fields[2].p_dof = fields[ej + 1][ei + 1].p_dof;
988: nodal_fields[3].u_dof = fields[ej][ei + 1].u_dof;
989: nodal_fields[3].v_dof = fields[ej][ei + 1].v_dof;
990: nodal_fields[3].p_dof = fields[ej][ei + 1].p_dof;
991: PetscFunctionReturn(PETSC_SUCCESS);
992: }
994: static PetscErrorCode DMDAIntegrateErrors(DM stokes_da, Vec X, Vec X_analytic)
995: {
996: DM cda;
997: Vec coords, X_analytic_local, X_local;
998: DMDACoor2d **_coords;
999: PetscInt sex, sey, mx, my;
1000: PetscInt ei, ej;
1001: PetscScalar el_coords[NODES_PER_EL * NSD];
1002: StokesDOF **stokes_analytic, **stokes;
1003: StokesDOF stokes_analytic_e[4], stokes_e[4];
1005: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
1006: PetscScalar Ni_p[NODES_PER_EL];
1007: PetscInt ngp;
1008: PetscScalar gp_xi[GAUSS_POINTS][2];
1009: PetscScalar gp_weight[GAUSS_POINTS];
1010: PetscInt p, i;
1011: PetscScalar J_p, fac;
1012: PetscScalar h, p_e_L2, u_e_L2, u_e_H1, p_L2, u_L2, u_H1, tp_L2, tu_L2, tu_H1;
1013: PetscInt M;
1014: PetscReal xymin[2], xymax[2];
1016: PetscFunctionBeginUser;
1017: /* define quadrature rule */
1018: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
1020: /* setup for coords */
1021: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
1022: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
1023: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1025: /* setup for analytic */
1026: PetscCall(DMCreateLocalVector(stokes_da, &X_analytic_local));
1027: PetscCall(DMGlobalToLocalBegin(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
1028: PetscCall(DMGlobalToLocalEnd(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
1029: PetscCall(DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic));
1031: /* setup for solution */
1032: PetscCall(DMCreateLocalVector(stokes_da, &X_local));
1033: PetscCall(DMGlobalToLocalBegin(stokes_da, X, INSERT_VALUES, X_local));
1034: PetscCall(DMGlobalToLocalEnd(stokes_da, X, INSERT_VALUES, X_local));
1035: PetscCall(DMDAVecGetArray(stokes_da, X_local, &stokes));
1037: PetscCall(DMDAGetInfo(stokes_da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1038: PetscCall(DMGetBoundingBox(stokes_da, xymin, xymax));
1040: h = (xymax[0] - xymin[0]) / ((PetscReal)M);
1042: tp_L2 = tu_L2 = tu_H1 = 0.0;
1044: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
1045: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
1046: for (ej = sey; ej < sey + my; ej++) {
1047: for (ei = sex; ei < sex + mx; ei++) {
1048: /* get coords for the element */
1049: PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
1050: PetscCall(StokesDAGetNodalFields(stokes, ei, ej, stokes_e));
1051: PetscCall(StokesDAGetNodalFields(stokes_analytic, ei, ej, stokes_analytic_e));
1053: /* evaluate integral */
1054: p_e_L2 = 0.0;
1055: u_e_L2 = 0.0;
1056: u_e_H1 = 0.0;
1057: for (p = 0; p < ngp; p++) {
1058: ConstructQ12D_Ni(gp_xi[p], Ni_p);
1059: ConstructQ12D_GNi(gp_xi[p], GNi_p);
1060: ConstructQ12D_GNx(GNi_p, GNx_p, el_coords, &J_p);
1061: fac = gp_weight[p] * J_p;
1063: for (i = 0; i < NODES_PER_EL; i++) {
1064: PetscScalar u_error, v_error;
1066: 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);
1068: u_error = stokes_e[i].u_dof - stokes_analytic_e[i].u_dof;
1069: v_error = stokes_e[i].v_dof - stokes_analytic_e[i].v_dof;
1070: u_e_L2 += fac * Ni_p[i] * (u_error * u_error + v_error * v_error);
1072: u_e_H1 = u_e_H1 + fac * (GNx_p[0][i] * u_error * GNx_p[0][i] * u_error /* du/dx */
1073: + GNx_p[1][i] * u_error * GNx_p[1][i] * u_error /* du/dy */
1074: + GNx_p[0][i] * v_error * GNx_p[0][i] * v_error /* dv/dx */
1075: + GNx_p[1][i] * v_error * GNx_p[1][i] * v_error); /* dv/dy */
1076: }
1077: }
1079: tp_L2 += p_e_L2;
1080: tu_L2 += u_e_L2;
1081: tu_H1 += u_e_H1;
1082: }
1083: }
1084: PetscCall(MPIU_Allreduce(&tp_L2, &p_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1085: PetscCall(MPIU_Allreduce(&tu_L2, &u_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1086: PetscCall(MPIU_Allreduce(&tu_H1, &u_H1, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1087: p_L2 = PetscSqrtScalar(p_L2);
1088: u_L2 = PetscSqrtScalar(u_L2);
1089: u_H1 = PetscSqrtScalar(u_H1);
1091: PetscCall(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)));
1093: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1095: PetscCall(DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic));
1096: PetscCall(VecDestroy(&X_analytic_local));
1097: PetscCall(DMDAVecRestoreArray(stokes_da, X_local, &stokes));
1098: PetscCall(VecDestroy(&X_local));
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: }
1102: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx, PetscInt my)
1103: {
1104: DM da_Stokes, da_prop;
1105: PetscInt u_dof, p_dof, dof, stencil_width;
1106: Mat A, B;
1107: DM prop_cda, vel_cda;
1108: Vec prop_coords, vel_coords;
1109: PetscInt si, sj, nx, ny, i, j, p;
1110: Vec f, X;
1111: PetscInt prop_dof, prop_stencil_width;
1112: Vec properties, l_properties;
1113: PetscReal dx, dy;
1114: PetscInt M, N;
1115: DMDACoor2d **_prop_coords, **_vel_coords;
1116: GaussPointCoefficients **element_props;
1117: PetscInt its;
1118: KSP ksp_S;
1119: PetscInt coefficient_structure = 0;
1120: PetscInt cpu_x, cpu_y, *lx = NULL, *ly = NULL;
1121: PetscBool use_gp_coords = PETSC_FALSE, set, output_gnuplot = PETSC_FALSE, glvis = PETSC_FALSE, change = PETSC_FALSE;
1122: char filename[PETSC_MAX_PATH_LEN];
1124: PetscFunctionBeginUser;
1126: PetscCall(PetscOptionsGetBool(NULL, NULL, "-gnuplot", &output_gnuplot, NULL));
1127: PetscCall(PetscOptionsGetBool(NULL, NULL, "-glvis", &glvis, NULL));
1128: PetscCall(PetscOptionsGetBool(NULL, NULL, "-change", &change, NULL));
1130: /* Generate the da for velocity and pressure */
1131: /*
1132: We use Q1 elements for the temperature.
1133: FEM has a 9-point stencil (BOX) or connectivity pattern
1134: Num nodes in each direction is mx+1, my+1
1135: */
1136: u_dof = U_DOFS; /* Vx, Vy - velocities */
1137: p_dof = P_DOFS; /* p - pressure */
1138: dof = u_dof + p_dof;
1139: stencil_width = 1;
1140: PetscCall(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));
1142: PetscCall(DMSetMatType(da_Stokes, MATAIJ));
1143: PetscCall(DMSetFromOptions(da_Stokes));
1144: PetscCall(DMSetUp(da_Stokes));
1145: PetscCall(DMDASetFieldName(da_Stokes, 0, "Vx"));
1146: PetscCall(DMDASetFieldName(da_Stokes, 1, "Vy"));
1147: PetscCall(DMDASetFieldName(da_Stokes, 2, "P"));
1149: /* unit box [0,1] x [0,1] */
1150: PetscCall(DMDASetUniformCoordinates(da_Stokes, 0.0, 1.0, 0.0, 1.0, 0., 0.));
1152: /* Generate element properties, we will assume all material properties are constant over the element */
1153: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1154: PetscCall(DMDAGetInfo(da_Stokes, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
1155: PetscCall(DMDAGetElementOwnershipRanges2d(da_Stokes, &lx, &ly));
1157: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
1158: prop_stencil_width = 0;
1159: PetscCall(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));
1160: PetscCall(DMSetFromOptions(da_prop));
1161: PetscCall(DMSetUp(da_prop));
1162: PetscCall(PetscFree(lx));
1163: PetscCall(PetscFree(ly));
1165: /* define centroid positions */
1166: PetscCall(DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1167: dx = 1.0 / ((PetscReal)(M));
1168: dy = 1.0 / ((PetscReal)(N));
1170: PetscCall(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));
1172: /* define coefficients */
1173: PetscCall(PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL));
1175: PetscCall(DMCreateGlobalVector(da_prop, &properties));
1176: PetscCall(DMCreateLocalVector(da_prop, &l_properties));
1177: PetscCall(DMDAVecGetArray(da_prop, l_properties, &element_props));
1179: PetscCall(DMGetCoordinateDM(da_prop, &prop_cda));
1180: PetscCall(DMGetCoordinatesLocal(da_prop, &prop_coords));
1181: PetscCall(DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords));
1183: PetscCall(DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0));
1185: PetscCall(DMGetCoordinateDM(da_Stokes, &vel_cda));
1186: PetscCall(DMGetCoordinatesLocal(da_Stokes, &vel_coords));
1187: PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));
1189: /* interpolate the coordinates */
1190: for (j = sj; j < sj + ny; j++) {
1191: for (i = si; i < si + nx; i++) {
1192: PetscInt ngp;
1193: PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
1194: PetscScalar el_coords[8];
1196: PetscCall(GetElementCoords(_vel_coords, i, j, el_coords));
1197: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
1199: for (p = 0; p < GAUSS_POINTS; p++) {
1200: PetscScalar gp_x, gp_y;
1201: PetscInt n;
1202: PetscScalar xi_p[2], Ni_p[4];
1204: xi_p[0] = gp_xi[p][0];
1205: xi_p[1] = gp_xi[p][1];
1206: ConstructQ12D_Ni(xi_p, Ni_p);
1208: gp_x = 0.0;
1209: gp_y = 0.0;
1210: for (n = 0; n < NODES_PER_EL; n++) {
1211: gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
1212: gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
1213: }
1214: element_props[j][i].gp_coords[2 * p] = gp_x;
1215: element_props[j][i].gp_coords[2 * p + 1] = gp_y;
1216: }
1217: }
1218: }
1220: /* define the coefficients */
1221: PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, NULL));
1223: for (j = sj; j < sj + ny; j++) {
1224: for (i = si; i < si + nx; i++) {
1225: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1226: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1227: PetscReal coord_x, coord_y;
1229: if (coefficient_structure == 0) {
1230: PetscReal opts_eta0, opts_eta1, opts_xc;
1231: PetscInt opts_nz;
1233: opts_eta0 = 1.0;
1234: opts_eta1 = 1.0;
1235: opts_xc = 0.5;
1236: opts_nz = 1;
1238: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL));
1239: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL));
1240: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL));
1241: PetscCall(PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL));
1243: for (p = 0; p < GAUSS_POINTS; p++) {
1244: coord_x = centroid_x;
1245: coord_y = centroid_y;
1246: if (use_gp_coords) {
1247: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1248: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1249: }
1251: element_props[j][i].eta[p] = opts_eta0;
1252: if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;
1254: element_props[j][i].fx[p] = 0.0;
1255: element_props[j][i].fy[p] = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1256: }
1257: } else if (coefficient_structure == 1) { /* square sinker */
1258: PetscReal opts_eta0, opts_eta1, opts_dx, opts_dy;
1260: opts_eta0 = 1.0;
1261: opts_eta1 = 1.0;
1262: opts_dx = 0.50;
1263: opts_dy = 0.50;
1265: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1266: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1267: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL));
1268: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL));
1270: for (p = 0; p < GAUSS_POINTS; p++) {
1271: coord_x = centroid_x;
1272: coord_y = centroid_y;
1273: if (use_gp_coords) {
1274: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1275: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1276: }
1278: element_props[j][i].eta[p] = opts_eta0;
1279: element_props[j][i].fx[p] = 0.0;
1280: element_props[j][i].fy[p] = 0.0;
1282: if ((coord_x > -0.5 * opts_dx + 0.5) && (coord_x < 0.5 * opts_dx + 0.5)) {
1283: if ((coord_y > -0.5 * opts_dy + 0.5) && (coord_y < 0.5 * opts_dy + 0.5)) {
1284: element_props[j][i].eta[p] = opts_eta1;
1285: element_props[j][i].fx[p] = 0.0;
1286: element_props[j][i].fy[p] = -1.0;
1287: }
1288: }
1289: }
1290: } else if (coefficient_structure == 2) { /* circular sinker */
1291: PetscReal opts_eta0, opts_eta1, opts_r, radius2;
1293: opts_eta0 = 1.0;
1294: opts_eta1 = 1.0;
1295: opts_r = 0.25;
1297: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1298: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1299: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL));
1301: for (p = 0; p < GAUSS_POINTS; p++) {
1302: coord_x = centroid_x;
1303: coord_y = centroid_y;
1304: if (use_gp_coords) {
1305: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1306: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1307: }
1309: element_props[j][i].eta[p] = opts_eta0;
1310: element_props[j][i].fx[p] = 0.0;
1311: element_props[j][i].fy[p] = 0.0;
1313: radius2 = (coord_x - 0.5) * (coord_x - 0.5) + (coord_y - 0.5) * (coord_y - 0.5);
1314: if (radius2 < opts_r * opts_r) {
1315: element_props[j][i].eta[p] = opts_eta1;
1316: element_props[j][i].fx[p] = 0.0;
1317: element_props[j][i].fy[p] = -1.0;
1318: }
1319: }
1320: } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1321: PetscReal opts_eta0, opts_eta1, opts_r, opts_dx, opts_dy, opts_c0x, opts_c0y, opts_s0x, opts_s0y, opts_phi, radius2;
1323: opts_eta0 = 1.0;
1324: opts_eta1 = 1.0;
1325: opts_r = 0.25;
1326: opts_c0x = 0.35; /* circle center */
1327: opts_c0y = 0.35;
1328: opts_s0x = 0.7; /* square center */
1329: opts_s0y = 0.7;
1330: opts_dx = 0.25;
1331: opts_dy = 0.25;
1332: opts_phi = 25;
1334: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1335: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1336: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL));
1337: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_c0x", &opts_c0x, NULL));
1338: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_c0y", &opts_c0y, NULL));
1339: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_s0x", &opts_s0x, NULL));
1340: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_s0y", &opts_s0y, NULL));
1341: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL));
1342: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL));
1343: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_phi", &opts_phi, NULL));
1344: opts_phi *= PETSC_PI / 180;
1346: for (p = 0; p < GAUSS_POINTS; p++) {
1347: coord_x = centroid_x;
1348: coord_y = centroid_y;
1349: if (use_gp_coords) {
1350: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1351: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1352: }
1354: element_props[j][i].eta[p] = opts_eta0;
1355: element_props[j][i].fx[p] = 0.0;
1356: element_props[j][i].fy[p] = -0.2;
1358: radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1359: 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)) {
1360: element_props[j][i].eta[p] = opts_eta1;
1361: element_props[j][i].fx[p] = 0.0;
1362: element_props[j][i].fy[p] = -1.0;
1363: }
1364: }
1365: } else if (coefficient_structure == 4) { /* subdomain jump */
1366: PetscReal opts_mag, opts_eta0;
1367: PetscInt opts_nz, px, py;
1368: PetscBool jump;
1370: opts_mag = 1.0;
1371: opts_eta0 = 1.0;
1372: opts_nz = 1;
1374: PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_eta0", &opts_eta0, NULL));
1375: PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_magnitude", &opts_mag, NULL));
1376: PetscCall(PetscOptionsGetInt(NULL, NULL, "-jump_nz", &opts_nz, NULL));
1377: PetscCall(DMDAGetInfo(da_Stokes, NULL, NULL, NULL, NULL, &px, &py, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1378: if (px % 2) {
1379: jump = (PetscBool)(PetscGlobalRank % 2);
1380: } else {
1381: jump = (PetscBool)((PetscGlobalRank / px) % 2 ? PetscGlobalRank % 2 : !(PetscGlobalRank % 2));
1382: }
1383: for (p = 0; p < GAUSS_POINTS; p++) {
1384: coord_x = centroid_x;
1385: coord_y = centroid_y;
1386: if (use_gp_coords) {
1387: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1388: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1389: }
1391: element_props[j][i].eta[p] = jump ? PetscPowReal(10.0, opts_mag) : opts_eta0;
1392: element_props[j][i].fx[p] = 0.0;
1393: element_props[j][i].fy[p] = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1394: }
1395: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
1396: }
1397: }
1398: PetscCall(DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords));
1400: PetscCall(DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords));
1402: PetscCall(DMDAVecRestoreArray(da_prop, l_properties, &element_props));
1403: PetscCall(DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties));
1404: PetscCall(DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties));
1406: if (output_gnuplot) {
1407: PetscCall(DMDACoordViewGnuplot2d(da_Stokes, "mesh"));
1408: PetscCall(DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for Stokes eqn.", "properties"));
1409: }
1411: if (glvis) {
1412: Vec glv_prop, etaf;
1413: PetscViewer view;
1414: PetscInt dim;
1415: const char *fec = {"FiniteElementCollection: L2_2D_P1"};
1417: PetscCall(DMGetDimension(da_Stokes, &dim));
1418: PetscCall(VecCreateSeq(PETSC_COMM_SELF, GAUSS_POINTS * mx * mx, &etaf));
1419: PetscCall(PetscObjectSetName((PetscObject)etaf, "viscosity"));
1420: PetscCall(PetscViewerGLVisOpen(PETSC_COMM_WORLD, PETSC_VIEWER_GLVIS_SOCKET, NULL, PETSC_DECIDE, &view));
1421: PetscCall(PetscViewerGLVisSetFields(view, 1, &fec, &dim, glvis_extract_eta, (PetscObject *)&etaf, da_prop, NULL));
1422: PetscCall(DMGetLocalVector(da_prop, &glv_prop));
1423: PetscCall(DMGlobalToLocalBegin(da_prop, properties, INSERT_VALUES, glv_prop));
1424: PetscCall(DMGlobalToLocalEnd(da_prop, properties, INSERT_VALUES, glv_prop));
1425: PetscCall(VecSetDM(etaf, da_Stokes));
1426: PetscCall(VecView(glv_prop, view));
1427: PetscCall(PetscViewerDestroy(&view));
1428: PetscCall(DMRestoreLocalVector(da_prop, &glv_prop));
1429: PetscCall(VecDestroy(&etaf));
1430: }
1432: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1433: PetscCall(DMCreateMatrix(da_Stokes, &A));
1434: PetscCall(DMCreateMatrix(da_Stokes, &B));
1435: PetscCall(DMCreateGlobalVector(da_Stokes, &f));
1436: PetscCall(DMCreateGlobalVector(da_Stokes, &X));
1438: /* assemble A11 */
1439: PetscCall(MatZeroEntries(A));
1440: PetscCall(MatZeroEntries(B));
1441: PetscCall(VecZeroEntries(f));
1443: PetscCall(AssembleA_Stokes(A, da_Stokes, da_prop, properties));
1444: PetscCall(AssembleA_PCStokes(B, da_Stokes, da_prop, properties));
1445: /* build force vector */
1446: PetscCall(AssembleF_Stokes(f, da_Stokes, da_prop, properties));
1448: PetscCall(DMDABCApplyFreeSlip(da_Stokes, A, f));
1449: PetscCall(DMDABCApplyFreeSlip(da_Stokes, B, NULL));
1451: /* SOLVE */
1452: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_S));
1453: PetscCall(KSPSetOptionsPrefix(ksp_S, "stokes_"));
1454: PetscCall(KSPSetDM(ksp_S, da_Stokes));
1455: PetscCall(KSPSetDMActive(ksp_S, PETSC_FALSE));
1456: PetscCall(KSPSetOperators(ksp_S, A, B));
1457: PetscCall(KSPSetFromOptions(ksp_S));
1458: {
1459: PC pc;
1460: const PetscInt ufields[] = {0, 1}, pfields[1] = {2};
1462: PetscCall(KSPGetPC(ksp_S, &pc));
1463: PetscCall(PCFieldSplitSetBlockSize(pc, 3));
1464: PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
1465: PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
1466: }
1468: {
1469: PC pc;
1470: PetscBool same = PETSC_FALSE;
1471: PetscCall(KSPGetPC(ksp_S, &pc));
1472: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &same));
1473: if (same) {
1474: PetscBool usedivmat = PETSC_FALSE;
1475: PetscCall(KSPSetOperators(ksp_S, A, A));
1477: PetscCall(PetscOptionsGetBool(NULL, NULL, "-stokes_pc_bddc_use_divergence", &usedivmat, NULL));
1478: if (usedivmat) {
1479: IS *fields, vel;
1480: PetscInt i, nf;
1482: PetscCall(DMCreateFieldDecomposition(da_Stokes, &nf, NULL, &fields, NULL));
1483: PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, fields, &vel));
1484: PetscCall(MatZeroRowsIS(B, fields[2], 1.0, NULL, NULL)); /* we put 1.0 on the diagonal to pick the pressure average too */
1485: PetscCall(MatTranspose(B, MAT_INPLACE_MATRIX, &B));
1486: PetscCall(MatZeroRowsIS(B, vel, 0.0, NULL, NULL));
1487: PetscCall(ISDestroy(&vel));
1488: PetscCall(PCBDDCSetDivergenceMat(pc, B, PETSC_FALSE, NULL));
1489: for (i = 0; i < nf; i++) PetscCall(ISDestroy(&fields[i]));
1490: PetscCall(PetscFree(fields));
1491: }
1492: }
1493: }
1495: {
1496: PC pc_S;
1497: KSP *sub_ksp, ksp_U;
1498: PetscInt nsplits;
1499: DM da_U;
1500: PetscBool is_pcfs;
1502: PetscCall(KSPSetUp(ksp_S));
1503: PetscCall(KSPGetPC(ksp_S, &pc_S));
1505: is_pcfs = PETSC_FALSE;
1506: PetscCall(PetscObjectTypeCompare((PetscObject)pc_S, PCFIELDSPLIT, &is_pcfs));
1508: if (is_pcfs) {
1509: PetscCall(PCFieldSplitGetSubKSP(pc_S, &nsplits, &sub_ksp));
1510: ksp_U = sub_ksp[0];
1511: PetscCall(PetscFree(sub_ksp));
1513: if (nsplits == 2) {
1514: PetscCall(DMDACreateCompatibleDMDA(da_Stokes, 2, &da_U));
1516: PetscCall(KSPSetDM(ksp_U, da_U));
1517: PetscCall(KSPSetDMActive(ksp_U, PETSC_FALSE));
1518: PetscCall(DMDestroy(&da_U));
1519: if (change) {
1520: const char opt[] = "-stokes_fieldsplit_u_pc_factor_mat_solver_type mumps";
1521: PetscCall(PetscOptionsInsertString(NULL, opt));
1522: PetscCall(KSPSetFromOptions(ksp_U));
1523: }
1524: PetscCall(KSPSetUp(ksp_U));
1525: }
1526: }
1527: }
1529: PetscCall(KSPSolve(ksp_S, f, X));
1531: PetscCall(PetscOptionsGetString(NULL, NULL, "-o", filename, sizeof(filename), &set));
1532: if (set) {
1533: char *ext;
1534: PetscViewer viewer;
1535: PetscBool flg;
1536: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1537: PetscCall(PetscStrrchr(filename, '.', &ext));
1538: PetscCall(PetscStrcmp("vts", ext, &flg));
1539: if (flg) {
1540: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1541: } else {
1542: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
1543: }
1544: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1545: PetscCall(PetscViewerFileSetName(viewer, filename));
1546: PetscCall(VecView(X, viewer));
1547: PetscCall(PetscViewerDestroy(&viewer));
1548: }
1549: if (output_gnuplot) PetscCall(DMDAViewGnuplot2d(da_Stokes, X, "Velocity solution for Stokes eqn.", "X"));
1551: if (glvis) {
1552: PetscViewer view;
1554: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &view));
1555: PetscCall(PetscViewerSetType(view, PETSCVIEWERGLVIS));
1556: PetscCall(VecView(X, view));
1557: PetscCall(PetscViewerDestroy(&view));
1558: }
1560: PetscCall(KSPGetIterationNumber(ksp_S, &its));
1562: if (coefficient_structure == 0) {
1563: PetscReal opts_eta0, opts_eta1, opts_xc;
1564: PetscInt opts_nz, N;
1565: DM da_Stokes_analytic;
1566: Vec X_analytic;
1567: PetscReal nrm1[3], nrm2[3], nrmI[3];
1569: opts_eta0 = 1.0;
1570: opts_eta1 = 1.0;
1571: opts_xc = 0.5;
1572: opts_nz = 1;
1574: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL));
1575: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL));
1576: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL));
1577: PetscCall(PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL));
1579: PetscCall(DMDACreateSolCx(opts_eta0, opts_eta1, opts_xc, opts_nz, mx, my, &da_Stokes_analytic, &X_analytic));
1580: if (output_gnuplot) PetscCall(DMDAViewGnuplot2d(da_Stokes_analytic, X_analytic, "Analytic solution for Stokes eqn.", "X_analytic"));
1581: PetscCall(DMDAIntegrateErrors(da_Stokes_analytic, X, X_analytic));
1583: PetscCall(VecAXPY(X_analytic, -1.0, X));
1584: PetscCall(VecGetSize(X_analytic, &N));
1585: N = N / 3;
1587: PetscCall(VecStrideNorm(X_analytic, 0, NORM_1, &nrm1[0]));
1588: PetscCall(VecStrideNorm(X_analytic, 0, NORM_2, &nrm2[0]));
1589: PetscCall(VecStrideNorm(X_analytic, 0, NORM_INFINITY, &nrmI[0]));
1591: PetscCall(VecStrideNorm(X_analytic, 1, NORM_1, &nrm1[1]));
1592: PetscCall(VecStrideNorm(X_analytic, 1, NORM_2, &nrm2[1]));
1593: PetscCall(VecStrideNorm(X_analytic, 1, NORM_INFINITY, &nrmI[1]));
1595: PetscCall(VecStrideNorm(X_analytic, 2, NORM_1, &nrm1[2]));
1596: PetscCall(VecStrideNorm(X_analytic, 2, NORM_2, &nrm2[2]));
1597: PetscCall(VecStrideNorm(X_analytic, 2, NORM_INFINITY, &nrmI[2]));
1599: PetscCall(DMDestroy(&da_Stokes_analytic));
1600: PetscCall(VecDestroy(&X_analytic));
1601: }
1603: PetscCall(KSPDestroy(&ksp_S));
1604: PetscCall(VecDestroy(&X));
1605: PetscCall(VecDestroy(&f));
1606: PetscCall(MatDestroy(&A));
1607: PetscCall(MatDestroy(&B));
1609: PetscCall(DMDestroy(&da_Stokes));
1610: PetscCall(DMDestroy(&da_prop));
1612: PetscCall(VecDestroy(&properties));
1613: PetscCall(VecDestroy(&l_properties));
1614: PetscFunctionReturn(PETSC_SUCCESS);
1615: }
1617: int main(int argc, char **args)
1618: {
1619: PetscInt mx, my;
1621: PetscFunctionBeginUser;
1622: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
1623: mx = my = 10;
1624: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1625: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1626: PetscCall(solve_stokes_2d_coupled(mx, my));
1627: PetscCall(PetscFinalize());
1628: return 0;
1629: }
1631: /* -------------------------- helpers for boundary conditions -------------------------------- */
1632: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
1633: {
1634: DM cda;
1635: Vec coords;
1636: PetscInt si, sj, nx, ny, i, j;
1637: PetscInt M, N;
1638: DMDACoor2d **_coords;
1639: const PetscInt *g_idx;
1640: PetscInt *bc_global_ids;
1641: PetscScalar *bc_vals;
1642: PetscInt nbcs;
1643: PetscInt n_dofs;
1644: ISLocalToGlobalMapping ltogm;
1646: PetscFunctionBeginUser;
1647: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1648: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1650: PetscCall(DMGetCoordinateDM(da, &cda));
1651: PetscCall(DMGetCoordinatesLocal(da, &coords));
1652: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1653: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1654: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1656: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1657: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1659: /* init the entries to -1 so VecSetValues will ignore them */
1660: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1662: i = nx - 1;
1663: for (j = 0; j < ny; j++) {
1664: PetscInt local_id;
1666: local_id = i + j * nx;
1668: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1670: bc_vals[j] = 0.0;
1671: }
1672: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1673: nbcs = 0;
1674: if ((si + nx) == (M)) nbcs = ny;
1676: if (b) {
1677: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1678: PetscCall(VecAssemblyBegin(b));
1679: PetscCall(VecAssemblyEnd(b));
1680: }
1681: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1683: PetscCall(PetscFree(bc_vals));
1684: PetscCall(PetscFree(bc_global_ids));
1686: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1687: PetscFunctionReturn(PETSC_SUCCESS);
1688: }
1690: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
1691: {
1692: DM cda;
1693: Vec coords;
1694: PetscInt si, sj, nx, ny, i, j;
1695: PetscInt M, N;
1696: DMDACoor2d **_coords;
1697: const PetscInt *g_idx;
1698: PetscInt *bc_global_ids;
1699: PetscScalar *bc_vals;
1700: PetscInt nbcs;
1701: PetscInt n_dofs;
1702: ISLocalToGlobalMapping ltogm;
1704: PetscFunctionBeginUser;
1705: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1706: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1708: PetscCall(DMGetCoordinateDM(da, &cda));
1709: PetscCall(DMGetCoordinatesLocal(da, &coords));
1710: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1711: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1712: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1714: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1715: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1717: /* init the entries to -1 so VecSetValues will ignore them */
1718: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1720: i = 0;
1721: for (j = 0; j < ny; j++) {
1722: PetscInt local_id;
1724: local_id = i + j * nx;
1726: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1728: bc_vals[j] = 0.0;
1729: }
1730: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1731: nbcs = 0;
1732: if (si == 0) nbcs = ny;
1734: if (b) {
1735: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1736: PetscCall(VecAssemblyBegin(b));
1737: PetscCall(VecAssemblyEnd(b));
1738: }
1740: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1742: PetscCall(PetscFree(bc_vals));
1743: PetscCall(PetscFree(bc_global_ids));
1745: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1746: PetscFunctionReturn(PETSC_SUCCESS);
1747: }
1749: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
1750: {
1751: DM cda;
1752: Vec coords;
1753: PetscInt si, sj, nx, ny, i, j;
1754: PetscInt M, N;
1755: DMDACoor2d **_coords;
1756: const PetscInt *g_idx;
1757: PetscInt *bc_global_ids;
1758: PetscScalar *bc_vals;
1759: PetscInt nbcs;
1760: PetscInt n_dofs;
1761: ISLocalToGlobalMapping ltogm;
1763: PetscFunctionBeginUser;
1764: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1765: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1767: PetscCall(DMGetCoordinateDM(da, &cda));
1768: PetscCall(DMGetCoordinatesLocal(da, &coords));
1769: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1770: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1771: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1773: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1774: PetscCall(PetscMalloc1(nx, &bc_vals));
1776: /* init the entries to -1 so VecSetValues will ignore them */
1777: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1779: j = ny - 1;
1780: for (i = 0; i < nx; i++) {
1781: PetscInt local_id;
1783: local_id = i + j * nx;
1785: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1787: bc_vals[i] = 0.0;
1788: }
1789: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1790: nbcs = 0;
1791: if ((sj + ny) == (N)) nbcs = nx;
1793: if (b) {
1794: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1795: PetscCall(VecAssemblyBegin(b));
1796: PetscCall(VecAssemblyEnd(b));
1797: }
1798: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1800: PetscCall(PetscFree(bc_vals));
1801: PetscCall(PetscFree(bc_global_ids));
1803: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1804: PetscFunctionReturn(PETSC_SUCCESS);
1805: }
1807: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
1808: {
1809: DM cda;
1810: Vec coords;
1811: PetscInt si, sj, nx, ny, i, j;
1812: PetscInt M, N;
1813: DMDACoor2d **_coords;
1814: const PetscInt *g_idx;
1815: PetscInt *bc_global_ids;
1816: PetscScalar *bc_vals;
1817: PetscInt nbcs;
1818: PetscInt n_dofs;
1819: ISLocalToGlobalMapping ltogm;
1821: PetscFunctionBeginUser;
1822: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1823: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1825: PetscCall(DMGetCoordinateDM(da, &cda));
1826: PetscCall(DMGetCoordinatesLocal(da, &coords));
1827: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1828: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1829: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1831: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1832: PetscCall(PetscMalloc1(nx, &bc_vals));
1834: /* init the entries to -1 so VecSetValues will ignore them */
1835: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1837: j = 0;
1838: for (i = 0; i < nx; i++) {
1839: PetscInt local_id;
1841: local_id = i + j * nx;
1843: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1845: bc_vals[i] = 0.0;
1846: }
1847: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1848: nbcs = 0;
1849: if (sj == 0) nbcs = nx;
1851: if (b) {
1852: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1853: PetscCall(VecAssemblyBegin(b));
1854: PetscCall(VecAssemblyEnd(b));
1855: }
1856: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1858: PetscCall(PetscFree(bc_vals));
1859: PetscCall(PetscFree(bc_global_ids));
1861: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1862: PetscFunctionReturn(PETSC_SUCCESS);
1863: }
1865: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
1866: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes, Mat A, Vec f)
1867: {
1868: PetscFunctionBeginUser;
1869: PetscCall(BCApplyZero_NORTH(da_Stokes, 1, A, f));
1870: PetscCall(BCApplyZero_EAST(da_Stokes, 0, A, f));
1871: PetscCall(BCApplyZero_SOUTH(da_Stokes, 1, A, f));
1872: PetscCall(BCApplyZero_WEST(da_Stokes, 0, A, f));
1873: PetscFunctionReturn(PETSC_SUCCESS);
1874: }
1876: /*TEST
1878: build:
1879: requires: !complex !single
1881: test:
1882: 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
1884: testset:
1885: 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
1886: test:
1887: suffix: 2
1888: args:
1889: output_file: output/ex43_1.out
1890: test:
1891: requires: mumps
1892: suffix: 2_mumps
1893: args: -change true -stokes_ksp_view
1894: output_file: output/ex43_2_mumps.out
1895: # mumps INFO,INFOG,RINFO,RINFOG may vary on different archs, so keep just a stable one
1896: filter: grep -E -v "(INFOG\([^7]|INFO\(|\[0\])" | grep -v " total: nonzeros="
1898: test:
1899: suffix: 3
1900: nsize: 4
1901: 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
1903: test:
1904: suffix: 4
1905: nsize: 4
1906: 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
1908: test:
1909: suffix: 5
1910: nsize: 4
1911: 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
1913: test:
1914: suffix: 6
1915: nsize: 8
1916: 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
1918: test:
1919: suffix: bjacobi
1920: nsize: 4
1921: args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type aij -stokes_ksp_converged_reason
1923: test:
1924: suffix: bjacobi_baij
1925: nsize: 4
1926: args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type baij -stokes_ksp_converged_reason
1927: output_file: output/ex43_bjacobi.out
1929: test:
1930: suffix: nested_gmg
1931: nsize: 4
1932: 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
1934: test:
1935: suffix: fetidp
1936: nsize: 8
1937: 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
1939: test:
1940: suffix: fetidp_unsym
1941: nsize: 8
1942: 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
1944: test:
1945: suffix: bddc_stokes_deluxe
1946: nsize: 8
1947: 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
1949: test:
1950: suffix: bddc_stokes_subdomainjump_deluxe
1951: nsize: 9
1952: 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
1954: TEST*/