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[], PetscCtx 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;
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 (PetscInt 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 (PetscInt 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;
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 (PetscInt 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: PetscFunctionBeginUser;
844: for (PetscInt 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: PetscFunctionReturn(PETSC_SUCCESS);
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 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;
870: PetscFunctionBeginUser;
871: /* setup for coords */
872: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
873: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
874: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
876: /* setup for coefficients */
877: PetscCall(DMGetLocalVector(properties_da, &local_properties));
878: PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
879: PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
880: PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));
882: /* get access to the vector */
883: PetscCall(DMGetLocalVector(stokes_da, &local_F));
884: PetscCall(VecZeroEntries(local_F));
885: PetscCall(DMDAVecGetArray(stokes_da, local_F, &ff));
887: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
888: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
889: for (ej = sey; ej < sey + my; ej++) {
890: for (PetscInt ei = sex; ei < sex + mx; ei++) {
891: /* get coords for the element */
892: PetscCall(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: PetscCall(PetscMemzero(Fe, sizeof(Fe)));
900: PetscCall(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: PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
908: PetscCall(DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, p_eqn, Fe, He));
909: }
910: }
912: PetscCall(DMDAVecRestoreArray(stokes_da, local_F, &ff));
913: PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
914: PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
915: PetscCall(DMRestoreLocalVector(stokes_da, &local_F));
917: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
919: PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
920: PetscCall(DMRestoreLocalVector(properties_da, &local_properties));
921: PetscFunctionReturn(PETSC_SUCCESS);
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;
933: PetscFunctionBeginUser;
934: 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));
935: PetscCall(DMSetFromOptions(da));
936: PetscCall(DMSetUp(da));
937: PetscCall(DMDASetFieldName(da, 0, "anlytic_Vx"));
938: PetscCall(DMDASetFieldName(da, 1, "anlytic_Vy"));
939: PetscCall(DMDASetFieldName(da, 2, "analytic_P"));
941: PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0., 0.));
943: PetscCall(DMGetCoordinatesLocal(da, &coords));
944: PetscCall(DMGetCoordinateDM(da, &cda));
945: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
947: PetscCall(DMCreateGlobalVector(da, &X));
948: PetscCall(DMDAVecGetArray(da, X, &_stokes));
950: PetscCall(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: PetscCall(DMDAVecRestoreArray(da, X, &_stokes));
966: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
968: *_da = da;
969: *_X = X;
970: PetscFunctionReturn(PETSC_SUCCESS);
971: }
973: static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields, PetscInt ei, PetscInt ej, StokesDOF nodal_fields[])
974: {
975: PetscFunctionBeginUser;
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: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscScalar el_coords[NODES_PER_EL * NSD];
999: StokesDOF **stokes_analytic, **stokes;
1000: StokesDOF stokes_analytic_e[4], stokes_e[4];
1002: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
1003: PetscScalar Ni_p[NODES_PER_EL];
1004: PetscInt ngp;
1005: PetscScalar gp_xi[GAUSS_POINTS][2];
1006: PetscScalar gp_weight[GAUSS_POINTS];
1007: PetscInt p, i;
1008: PetscScalar J_p, fac;
1009: PetscScalar h, p_e_L2, u_e_L2, u_e_H1, p_L2, u_L2, u_H1, tp_L2, tu_L2, tu_H1;
1010: PetscInt M;
1011: PetscReal xymin[2], xymax[2];
1013: PetscFunctionBeginUser;
1014: /* define quadrature rule */
1015: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
1017: /* setup for coords */
1018: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
1019: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
1020: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1022: /* setup for analytic */
1023: PetscCall(DMCreateLocalVector(stokes_da, &X_analytic_local));
1024: PetscCall(DMGlobalToLocalBegin(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
1025: PetscCall(DMGlobalToLocalEnd(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
1026: PetscCall(DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic));
1028: /* setup for solution */
1029: PetscCall(DMCreateLocalVector(stokes_da, &X_local));
1030: PetscCall(DMGlobalToLocalBegin(stokes_da, X, INSERT_VALUES, X_local));
1031: PetscCall(DMGlobalToLocalEnd(stokes_da, X, INSERT_VALUES, X_local));
1032: PetscCall(DMDAVecGetArray(stokes_da, X_local, &stokes));
1034: PetscCall(DMDAGetInfo(stokes_da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1035: PetscCall(DMGetBoundingBox(stokes_da, xymin, xymax));
1037: h = (xymax[0] - xymin[0]) / ((PetscReal)M);
1039: tp_L2 = tu_L2 = tu_H1 = 0.0;
1041: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
1042: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
1043: for (PetscInt ej = sey; ej < sey + my; ej++) {
1044: for (PetscInt ei = sex; ei < sex + mx; ei++) {
1045: /* get coords for the element */
1046: PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
1047: PetscCall(StokesDAGetNodalFields(stokes, ei, ej, stokes_e));
1048: PetscCall(StokesDAGetNodalFields(stokes_analytic, ei, ej, stokes_analytic_e));
1050: /* evaluate integral */
1051: p_e_L2 = 0.0;
1052: u_e_L2 = 0.0;
1053: u_e_H1 = 0.0;
1054: for (p = 0; p < ngp; p++) {
1055: ConstructQ12D_Ni(gp_xi[p], Ni_p);
1056: ConstructQ12D_GNi(gp_xi[p], GNi_p);
1057: ConstructQ12D_GNx(GNi_p, GNx_p, el_coords, &J_p);
1058: fac = gp_weight[p] * J_p;
1060: for (i = 0; i < NODES_PER_EL; i++) {
1061: PetscScalar u_error, v_error;
1063: 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);
1065: u_error = stokes_e[i].u_dof - stokes_analytic_e[i].u_dof;
1066: v_error = stokes_e[i].v_dof - stokes_analytic_e[i].v_dof;
1067: u_e_L2 += fac * Ni_p[i] * (u_error * u_error + v_error * v_error);
1069: u_e_H1 = u_e_H1 + fac * (GNx_p[0][i] * u_error * GNx_p[0][i] * u_error /* du/dx */
1070: + GNx_p[1][i] * u_error * GNx_p[1][i] * u_error /* du/dy */
1071: + GNx_p[0][i] * v_error * GNx_p[0][i] * v_error /* dv/dx */
1072: + GNx_p[1][i] * v_error * GNx_p[1][i] * v_error); /* dv/dy */
1073: }
1074: }
1076: tp_L2 += p_e_L2;
1077: tu_L2 += u_e_L2;
1078: tu_H1 += u_e_H1;
1079: }
1080: }
1081: PetscCallMPI(MPIU_Allreduce(&tp_L2, &p_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1082: PetscCallMPI(MPIU_Allreduce(&tu_L2, &u_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1083: PetscCallMPI(MPIU_Allreduce(&tu_H1, &u_H1, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1084: p_L2 = PetscSqrtScalar(p_L2);
1085: u_L2 = PetscSqrtScalar(u_L2);
1086: u_H1 = PetscSqrtScalar(u_H1);
1088: 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)));
1090: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1092: PetscCall(DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic));
1093: PetscCall(VecDestroy(&X_analytic_local));
1094: PetscCall(DMDAVecRestoreArray(stokes_da, X_local, &stokes));
1095: PetscCall(VecDestroy(&X_local));
1096: PetscFunctionReturn(PETSC_SUCCESS);
1097: }
1099: static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx, PetscInt my)
1100: {
1101: DM da_Stokes, da_prop;
1102: PetscInt u_dof, p_dof, dof, stencil_width;
1103: Mat A, B;
1104: DM prop_cda, vel_cda;
1105: Vec prop_coords, vel_coords;
1106: PetscInt si, sj, nx, ny, i, j, p;
1107: Vec f, X;
1108: PetscInt prop_dof, prop_stencil_width;
1109: Vec properties, l_properties;
1110: PetscReal dx, dy;
1111: PetscInt M, N;
1112: DMDACoor2d **_prop_coords, **_vel_coords;
1113: GaussPointCoefficients **element_props;
1114: PetscInt its;
1115: KSP ksp_S;
1116: PetscInt coefficient_structure = 0;
1117: PetscInt cpu_x, cpu_y, *lx = NULL, *ly = NULL;
1118: PetscBool use_gp_coords = PETSC_FALSE, set, output_gnuplot = PETSC_FALSE, glvis = PETSC_FALSE, change = PETSC_FALSE;
1119: char filename[PETSC_MAX_PATH_LEN];
1121: PetscFunctionBeginUser;
1122: PetscCall(PetscOptionsGetBool(NULL, NULL, "-gnuplot", &output_gnuplot, NULL));
1123: PetscCall(PetscOptionsGetBool(NULL, NULL, "-glvis", &glvis, NULL));
1124: PetscCall(PetscOptionsGetBool(NULL, NULL, "-change", &change, NULL));
1126: /* Generate the da for velocity and pressure */
1127: /*
1128: We use Q1 elements for the temperature.
1129: FEM has a 9-point stencil (BOX) or connectivity pattern
1130: Num nodes in each direction is mx+1, my+1
1131: */
1132: u_dof = U_DOFS; /* Vx, Vy - velocities */
1133: p_dof = P_DOFS; /* p - pressure */
1134: dof = u_dof + p_dof;
1135: stencil_width = 1;
1136: 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));
1138: PetscCall(DMSetMatType(da_Stokes, MATAIJ));
1139: PetscCall(DMSetFromOptions(da_Stokes));
1140: PetscCall(DMSetUp(da_Stokes));
1141: PetscCall(DMDASetFieldName(da_Stokes, 0, "Vx"));
1142: PetscCall(DMDASetFieldName(da_Stokes, 1, "Vy"));
1143: PetscCall(DMDASetFieldName(da_Stokes, 2, "P"));
1145: /* unit box [0,1] x [0,1] */
1146: PetscCall(DMDASetUniformCoordinates(da_Stokes, 0.0, 1.0, 0.0, 1.0, 0., 0.));
1148: /* Generate element properties, we will assume all material properties are constant over the element */
1149: /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1150: PetscCall(DMDAGetInfo(da_Stokes, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
1151: PetscCall(DMDAGetElementOwnershipRanges2d(da_Stokes, &lx, &ly));
1153: prop_dof = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
1154: prop_stencil_width = 0;
1155: 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));
1156: PetscCall(DMSetFromOptions(da_prop));
1157: PetscCall(DMSetUp(da_prop));
1158: PetscCall(PetscFree(lx));
1159: PetscCall(PetscFree(ly));
1161: /* define centroid positions */
1162: PetscCall(DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1163: dx = 1.0 / (PetscReal)M;
1164: dy = 1.0 / (PetscReal)N;
1166: 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));
1168: /* define coefficients */
1169: PetscCall(PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL));
1171: PetscCall(DMCreateGlobalVector(da_prop, &properties));
1172: PetscCall(DMCreateLocalVector(da_prop, &l_properties));
1173: PetscCall(DMDAVecGetArray(da_prop, l_properties, &element_props));
1175: PetscCall(DMGetCoordinateDM(da_prop, &prop_cda));
1176: PetscCall(DMGetCoordinatesLocal(da_prop, &prop_coords));
1177: PetscCall(DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords));
1179: PetscCall(DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0));
1181: PetscCall(DMGetCoordinateDM(da_Stokes, &vel_cda));
1182: PetscCall(DMGetCoordinatesLocal(da_Stokes, &vel_coords));
1183: PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));
1185: /* interpolate the coordinates */
1186: for (j = sj; j < sj + ny; j++) {
1187: for (i = si; i < si + nx; i++) {
1188: PetscInt ngp;
1189: PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
1190: PetscScalar el_coords[8];
1192: PetscCall(GetElementCoords(_vel_coords, i, j, el_coords));
1193: ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
1195: for (p = 0; p < GAUSS_POINTS; p++) {
1196: PetscScalar gp_x, gp_y;
1197: PetscScalar xi_p[2], Ni_p[4];
1199: xi_p[0] = gp_xi[p][0];
1200: xi_p[1] = gp_xi[p][1];
1201: ConstructQ12D_Ni(xi_p, Ni_p);
1203: gp_x = 0.0;
1204: gp_y = 0.0;
1205: for (PetscInt n = 0; n < NODES_PER_EL; n++) {
1206: gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
1207: gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
1208: }
1209: element_props[j][i].gp_coords[2 * p] = gp_x;
1210: element_props[j][i].gp_coords[2 * p + 1] = gp_y;
1211: }
1212: }
1213: }
1215: /* define the coefficients */
1216: PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, NULL));
1218: for (j = sj; j < sj + ny; j++) {
1219: for (i = si; i < si + nx; i++) {
1220: PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1221: PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1222: PetscReal coord_x, coord_y;
1224: if (coefficient_structure == 0) {
1225: PetscReal opts_eta0, opts_eta1, opts_xc;
1226: PetscInt opts_nz;
1228: opts_eta0 = 1.0;
1229: opts_eta1 = 1.0;
1230: opts_xc = 0.5;
1231: opts_nz = 1;
1233: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL));
1234: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL));
1235: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL));
1236: PetscCall(PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL));
1238: for (p = 0; p < GAUSS_POINTS; p++) {
1239: coord_x = centroid_x;
1240: coord_y = centroid_y;
1241: if (use_gp_coords) {
1242: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1243: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1244: }
1246: element_props[j][i].eta[p] = opts_eta0;
1247: if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;
1249: element_props[j][i].fx[p] = 0.0;
1250: element_props[j][i].fy[p] = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1251: }
1252: } else if (coefficient_structure == 1) { /* square sinker */
1253: PetscReal opts_eta0, opts_eta1, opts_dx, opts_dy;
1255: opts_eta0 = 1.0;
1256: opts_eta1 = 1.0;
1257: opts_dx = 0.50;
1258: opts_dy = 0.50;
1260: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1261: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1262: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL));
1263: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL));
1265: for (p = 0; p < GAUSS_POINTS; p++) {
1266: coord_x = centroid_x;
1267: coord_y = centroid_y;
1268: if (use_gp_coords) {
1269: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1270: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1271: }
1273: element_props[j][i].eta[p] = opts_eta0;
1274: element_props[j][i].fx[p] = 0.0;
1275: element_props[j][i].fy[p] = 0.0;
1277: if ((coord_x > -0.5 * opts_dx + 0.5) && (coord_x < 0.5 * opts_dx + 0.5)) {
1278: if ((coord_y > -0.5 * opts_dy + 0.5) && (coord_y < 0.5 * opts_dy + 0.5)) {
1279: element_props[j][i].eta[p] = opts_eta1;
1280: element_props[j][i].fx[p] = 0.0;
1281: element_props[j][i].fy[p] = -1.0;
1282: }
1283: }
1284: }
1285: } else if (coefficient_structure == 2) { /* circular sinker */
1286: PetscReal opts_eta0, opts_eta1, opts_r, radius2;
1288: opts_eta0 = 1.0;
1289: opts_eta1 = 1.0;
1290: opts_r = 0.25;
1292: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1293: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1294: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL));
1296: for (p = 0; p < GAUSS_POINTS; p++) {
1297: coord_x = centroid_x;
1298: coord_y = centroid_y;
1299: if (use_gp_coords) {
1300: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1301: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1302: }
1304: element_props[j][i].eta[p] = opts_eta0;
1305: element_props[j][i].fx[p] = 0.0;
1306: element_props[j][i].fy[p] = 0.0;
1308: radius2 = (coord_x - 0.5) * (coord_x - 0.5) + (coord_y - 0.5) * (coord_y - 0.5);
1309: if (radius2 < opts_r * opts_r) {
1310: element_props[j][i].eta[p] = opts_eta1;
1311: element_props[j][i].fx[p] = 0.0;
1312: element_props[j][i].fy[p] = -1.0;
1313: }
1314: }
1315: } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1316: PetscReal opts_eta0, opts_eta1, opts_r, opts_dx, opts_dy, opts_c0x, opts_c0y, opts_s0x, opts_s0y, opts_phi, radius2;
1318: opts_eta0 = 1.0;
1319: opts_eta1 = 1.0;
1320: opts_r = 0.25;
1321: opts_c0x = 0.35; /* circle center */
1322: opts_c0y = 0.35;
1323: opts_s0x = 0.7; /* square center */
1324: opts_s0y = 0.7;
1325: opts_dx = 0.25;
1326: opts_dy = 0.25;
1327: opts_phi = 25;
1329: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1330: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1331: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL));
1332: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_c0x", &opts_c0x, NULL));
1333: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_c0y", &opts_c0y, NULL));
1334: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_s0x", &opts_s0x, NULL));
1335: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_s0y", &opts_s0y, NULL));
1336: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL));
1337: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL));
1338: PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_phi", &opts_phi, NULL));
1339: opts_phi *= PETSC_PI / 180;
1341: for (p = 0; p < GAUSS_POINTS; p++) {
1342: coord_x = centroid_x;
1343: coord_y = centroid_y;
1344: if (use_gp_coords) {
1345: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1346: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1347: }
1349: element_props[j][i].eta[p] = opts_eta0;
1350: element_props[j][i].fx[p] = 0.0;
1351: element_props[j][i].fy[p] = -0.2;
1353: radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1354: 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)) {
1355: element_props[j][i].eta[p] = opts_eta1;
1356: element_props[j][i].fx[p] = 0.0;
1357: element_props[j][i].fy[p] = -1.0;
1358: }
1359: }
1360: } else if (coefficient_structure == 4) { /* subdomain jump */
1361: PetscReal opts_mag, opts_eta0;
1362: PetscInt opts_nz, px, py;
1363: PetscBool jump;
1365: opts_mag = 1.0;
1366: opts_eta0 = 1.0;
1367: opts_nz = 1;
1369: PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_eta0", &opts_eta0, NULL));
1370: PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_magnitude", &opts_mag, NULL));
1371: PetscCall(PetscOptionsGetInt(NULL, NULL, "-jump_nz", &opts_nz, NULL));
1372: PetscCall(DMDAGetInfo(da_Stokes, NULL, NULL, NULL, NULL, &px, &py, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1373: if (px % 2) {
1374: jump = (PetscBool)(PetscGlobalRank % 2);
1375: } else {
1376: jump = (PetscBool)((PetscGlobalRank / px) % 2 ? PetscGlobalRank % 2 : !(PetscGlobalRank % 2));
1377: }
1378: for (p = 0; p < GAUSS_POINTS; p++) {
1379: coord_x = centroid_x;
1380: coord_y = centroid_y;
1381: if (use_gp_coords) {
1382: coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1383: coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1384: }
1386: element_props[j][i].eta[p] = jump ? PetscPowReal(10.0, opts_mag) : opts_eta0;
1387: element_props[j][i].fx[p] = 0.0;
1388: element_props[j][i].fy[p] = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1389: }
1390: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
1391: }
1392: }
1393: PetscCall(DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords));
1395: PetscCall(DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords));
1397: PetscCall(DMDAVecRestoreArray(da_prop, l_properties, &element_props));
1398: PetscCall(DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties));
1399: PetscCall(DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties));
1401: if (output_gnuplot) {
1402: PetscCall(DMDACoordViewGnuplot2d(da_Stokes, "mesh"));
1403: PetscCall(DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for Stokes eqn.", "properties"));
1404: }
1406: if (glvis) {
1407: Vec glv_prop, etaf;
1408: PetscViewer view;
1409: PetscInt dim;
1410: const char *fec = {"FiniteElementCollection: L2_2D_P1"};
1412: PetscCall(DMGetDimension(da_Stokes, &dim));
1413: PetscCall(VecCreateSeq(PETSC_COMM_SELF, GAUSS_POINTS * mx * mx, &etaf));
1414: PetscCall(PetscObjectSetName((PetscObject)etaf, "viscosity"));
1415: PetscCall(PetscViewerGLVisOpen(PETSC_COMM_WORLD, PETSC_VIEWER_GLVIS_SOCKET, NULL, PETSC_DECIDE, &view));
1416: PetscCall(PetscViewerGLVisSetFields(view, 1, &fec, &dim, glvis_extract_eta, (PetscObject *)&etaf, da_prop, NULL));
1417: PetscCall(DMGetLocalVector(da_prop, &glv_prop));
1418: PetscCall(DMGlobalToLocalBegin(da_prop, properties, INSERT_VALUES, glv_prop));
1419: PetscCall(DMGlobalToLocalEnd(da_prop, properties, INSERT_VALUES, glv_prop));
1420: PetscCall(VecSetDM(etaf, da_Stokes));
1421: PetscCall(VecView(glv_prop, view));
1422: PetscCall(PetscViewerDestroy(&view));
1423: PetscCall(DMRestoreLocalVector(da_prop, &glv_prop));
1424: PetscCall(VecDestroy(&etaf));
1425: }
1427: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1428: PetscCall(DMCreateMatrix(da_Stokes, &A));
1429: PetscCall(DMCreateMatrix(da_Stokes, &B));
1430: PetscCall(DMCreateGlobalVector(da_Stokes, &f));
1431: PetscCall(DMCreateGlobalVector(da_Stokes, &X));
1433: /* assemble A11 */
1434: PetscCall(MatZeroEntries(A));
1435: PetscCall(MatZeroEntries(B));
1437: PetscCall(AssembleA_Stokes(A, da_Stokes, da_prop, properties));
1438: PetscCall(AssembleA_PCStokes(B, da_Stokes, da_prop, properties));
1439: /* build force vector */
1440: PetscCall(AssembleF_Stokes(f, da_Stokes, da_prop, properties));
1442: PetscCall(DMDABCApplyFreeSlip(da_Stokes, A, f));
1443: PetscCall(DMDABCApplyFreeSlip(da_Stokes, B, NULL));
1445: /* SOLVE */
1446: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_S));
1447: PetscCall(KSPSetOptionsPrefix(ksp_S, "stokes_"));
1448: PetscCall(KSPSetDM(ksp_S, da_Stokes));
1449: PetscCall(KSPSetDMActive(ksp_S, KSP_DMACTIVE_ALL, PETSC_FALSE));
1450: PetscCall(KSPSetOperators(ksp_S, A, B));
1451: PetscCall(KSPSetFromOptions(ksp_S));
1452: {
1453: PC pc;
1454: const PetscInt ufields[] = {0, 1}, pfields[1] = {2};
1456: PetscCall(KSPGetPC(ksp_S, &pc));
1457: PetscCall(PCFieldSplitSetBlockSize(pc, 3));
1458: PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
1459: PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
1460: }
1462: {
1463: PC pc;
1464: PetscBool same = PETSC_FALSE;
1465: PetscCall(KSPGetPC(ksp_S, &pc));
1466: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &same));
1467: if (same) {
1468: PetscBool usedivmat = PETSC_FALSE;
1469: PetscCall(KSPSetOperators(ksp_S, A, A));
1471: PetscCall(PetscOptionsGetBool(NULL, NULL, "-stokes_pc_bddc_use_divergence", &usedivmat, NULL));
1472: if (usedivmat) {
1473: IS *fields, vel;
1474: PetscInt nf;
1476: PetscCall(DMCreateFieldDecomposition(da_Stokes, &nf, NULL, &fields, NULL));
1477: PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, fields, &vel));
1478: PetscCall(MatZeroRowsIS(B, fields[2], 1.0, NULL, NULL)); /* we put 1.0 on the diagonal to pick the pressure average too */
1479: PetscCall(MatTranspose(B, MAT_INPLACE_MATRIX, &B));
1480: PetscCall(MatZeroRowsIS(B, vel, 0.0, NULL, NULL));
1481: PetscCall(ISDestroy(&vel));
1482: PetscCall(PCBDDCSetDivergenceMat(pc, B, PETSC_FALSE, NULL));
1483: for (PetscInt i = 0; i < nf; i++) PetscCall(ISDestroy(&fields[i]));
1484: PetscCall(PetscFree(fields));
1485: }
1486: }
1487: }
1489: {
1490: PC pc_S;
1491: KSP *sub_ksp, ksp_U;
1492: PetscInt nsplits;
1493: DM da_U;
1494: PetscBool is_pcfs;
1496: PetscCall(KSPSetUp(ksp_S));
1497: PetscCall(KSPGetPC(ksp_S, &pc_S));
1499: is_pcfs = PETSC_FALSE;
1500: PetscCall(PetscObjectTypeCompare((PetscObject)pc_S, PCFIELDSPLIT, &is_pcfs));
1502: if (is_pcfs) {
1503: PetscCall(PCFieldSplitGetSubKSP(pc_S, &nsplits, &sub_ksp));
1504: ksp_U = sub_ksp[0];
1505: PetscCall(PetscFree(sub_ksp));
1507: if (nsplits == 2) {
1508: PetscCall(DMDACreateCompatibleDMDA(da_Stokes, 2, &da_U));
1510: PetscCall(KSPSetDM(ksp_U, da_U));
1511: PetscCall(KSPSetDMActive(ksp_U, KSP_DMACTIVE_ALL, PETSC_FALSE));
1512: PetscCall(DMDestroy(&da_U));
1513: if (change) {
1514: const char opt[] = "-stokes_fieldsplit_u_pc_factor_mat_solver_type mumps";
1515: PetscCall(PetscOptionsInsertString(NULL, opt));
1516: PetscCall(KSPSetFromOptions(ksp_U));
1517: }
1518: PetscCall(KSPSetUp(ksp_U));
1519: }
1520: }
1521: }
1523: PetscCall(KSPSolve(ksp_S, f, X));
1525: PetscCall(PetscOptionsGetString(NULL, NULL, "-o", filename, sizeof(filename), &set));
1526: if (set) {
1527: char *ext;
1528: PetscViewer viewer;
1529: PetscBool flg;
1530: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1531: PetscCall(PetscStrrchr(filename, '.', &ext));
1532: PetscCall(PetscStrcmp("vts", ext, &flg));
1533: if (flg) {
1534: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1535: } else {
1536: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
1537: }
1538: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1539: PetscCall(PetscViewerFileSetName(viewer, filename));
1540: PetscCall(VecView(X, viewer));
1541: PetscCall(PetscViewerDestroy(&viewer));
1542: }
1543: if (output_gnuplot) PetscCall(DMDAViewGnuplot2d(da_Stokes, X, "Velocity solution for Stokes eqn.", "X"));
1545: if (glvis) {
1546: PetscViewer view;
1548: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &view));
1549: PetscCall(PetscViewerSetType(view, PETSCVIEWERGLVIS));
1550: PetscCall(VecView(X, view));
1551: PetscCall(PetscViewerDestroy(&view));
1552: }
1554: PetscCall(KSPGetIterationNumber(ksp_S, &its));
1556: if (coefficient_structure == 0) {
1557: PetscReal opts_eta0, opts_eta1, opts_xc;
1558: PetscInt opts_nz, N;
1559: DM da_Stokes_analytic;
1560: Vec X_analytic;
1561: PetscReal nrm1[3], nrm2[3], nrmI[3];
1563: opts_eta0 = 1.0;
1564: opts_eta1 = 1.0;
1565: opts_xc = 0.5;
1566: opts_nz = 1;
1568: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL));
1569: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL));
1570: PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL));
1571: PetscCall(PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL));
1573: PetscCall(DMDACreateSolCx(opts_eta0, opts_eta1, opts_xc, opts_nz, mx, my, &da_Stokes_analytic, &X_analytic));
1574: if (output_gnuplot) PetscCall(DMDAViewGnuplot2d(da_Stokes_analytic, X_analytic, "Analytic solution for Stokes eqn.", "X_analytic"));
1575: PetscCall(DMDAIntegrateErrors(da_Stokes_analytic, X, X_analytic));
1577: PetscCall(VecAXPY(X_analytic, -1.0, X));
1578: PetscCall(VecGetSize(X_analytic, &N));
1579: N = N / 3;
1581: PetscCall(VecStrideNorm(X_analytic, 0, NORM_1, &nrm1[0]));
1582: PetscCall(VecStrideNorm(X_analytic, 0, NORM_2, &nrm2[0]));
1583: PetscCall(VecStrideNorm(X_analytic, 0, NORM_INFINITY, &nrmI[0]));
1585: PetscCall(VecStrideNorm(X_analytic, 1, NORM_1, &nrm1[1]));
1586: PetscCall(VecStrideNorm(X_analytic, 1, NORM_2, &nrm2[1]));
1587: PetscCall(VecStrideNorm(X_analytic, 1, NORM_INFINITY, &nrmI[1]));
1589: PetscCall(VecStrideNorm(X_analytic, 2, NORM_1, &nrm1[2]));
1590: PetscCall(VecStrideNorm(X_analytic, 2, NORM_2, &nrm2[2]));
1591: PetscCall(VecStrideNorm(X_analytic, 2, NORM_INFINITY, &nrmI[2]));
1593: PetscCall(DMDestroy(&da_Stokes_analytic));
1594: PetscCall(VecDestroy(&X_analytic));
1595: }
1597: PetscCall(KSPDestroy(&ksp_S));
1598: PetscCall(VecDestroy(&X));
1599: PetscCall(VecDestroy(&f));
1600: PetscCall(MatDestroy(&A));
1601: PetscCall(MatDestroy(&B));
1603: PetscCall(DMDestroy(&da_Stokes));
1604: PetscCall(DMDestroy(&da_prop));
1606: PetscCall(VecDestroy(&properties));
1607: PetscCall(VecDestroy(&l_properties));
1608: PetscFunctionReturn(PETSC_SUCCESS);
1609: }
1611: int main(int argc, char **args)
1612: {
1613: PetscInt mx, my;
1615: PetscFunctionBeginUser;
1616: PetscCall(PetscInitialize(&argc, &args, NULL, help));
1617: mx = my = 10;
1618: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1619: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1620: PetscCall(solve_stokes_2d_coupled(mx, my));
1621: PetscCall(PetscFinalize());
1622: return 0;
1623: }
1625: /* -------------------------- helpers for boundary conditions -------------------------------- */
1626: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
1627: {
1628: DM cda;
1629: Vec coords;
1630: PetscInt si, sj, nx, ny, i, j;
1631: PetscInt M, N;
1632: DMDACoor2d **_coords;
1633: const PetscInt *g_idx;
1634: PetscInt *bc_global_ids;
1635: PetscScalar *bc_vals;
1636: PetscInt nbcs;
1637: PetscInt n_dofs;
1638: ISLocalToGlobalMapping ltogm;
1640: PetscFunctionBeginUser;
1641: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1642: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1644: PetscCall(DMGetCoordinateDM(da, &cda));
1645: PetscCall(DMGetCoordinatesLocal(da, &coords));
1646: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1647: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1648: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1650: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1651: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1653: /* init the entries to -1 so VecSetValues will ignore them */
1654: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1656: i = nx - 1;
1657: for (j = 0; j < ny; j++) {
1658: PetscInt local_id;
1660: local_id = i + j * nx;
1662: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1664: bc_vals[j] = 0.0;
1665: }
1666: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1667: nbcs = 0;
1668: if ((si + nx) == (M)) nbcs = ny;
1670: if (b) {
1671: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1672: PetscCall(VecAssemblyBegin(b));
1673: PetscCall(VecAssemblyEnd(b));
1674: }
1675: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1677: PetscCall(PetscFree(bc_vals));
1678: PetscCall(PetscFree(bc_global_ids));
1680: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1681: PetscFunctionReturn(PETSC_SUCCESS);
1682: }
1684: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
1685: {
1686: DM cda;
1687: Vec coords;
1688: PetscInt si, sj, nx, ny, i, j;
1689: PetscInt M, N;
1690: DMDACoor2d **_coords;
1691: const PetscInt *g_idx;
1692: PetscInt *bc_global_ids;
1693: PetscScalar *bc_vals;
1694: PetscInt nbcs;
1695: PetscInt n_dofs;
1696: ISLocalToGlobalMapping ltogm;
1698: PetscFunctionBeginUser;
1699: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1700: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1702: PetscCall(DMGetCoordinateDM(da, &cda));
1703: PetscCall(DMGetCoordinatesLocal(da, &coords));
1704: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1705: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1706: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1708: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1709: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1711: /* init the entries to -1 so VecSetValues will ignore them */
1712: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1714: i = 0;
1715: for (j = 0; j < ny; j++) {
1716: PetscInt local_id;
1718: local_id = i + j * nx;
1720: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1722: bc_vals[j] = 0.0;
1723: }
1724: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1725: nbcs = 0;
1726: if (si == 0) nbcs = ny;
1728: if (b) {
1729: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1730: PetscCall(VecAssemblyBegin(b));
1731: PetscCall(VecAssemblyEnd(b));
1732: }
1734: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1736: PetscCall(PetscFree(bc_vals));
1737: PetscCall(PetscFree(bc_global_ids));
1739: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1740: PetscFunctionReturn(PETSC_SUCCESS);
1741: }
1743: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
1744: {
1745: DM cda;
1746: Vec coords;
1747: PetscInt si, sj, nx, ny, i, j;
1748: PetscInt M, N;
1749: DMDACoor2d **_coords;
1750: const PetscInt *g_idx;
1751: PetscInt *bc_global_ids;
1752: PetscScalar *bc_vals;
1753: PetscInt nbcs;
1754: PetscInt n_dofs;
1755: ISLocalToGlobalMapping ltogm;
1757: PetscFunctionBeginUser;
1758: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1759: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1761: PetscCall(DMGetCoordinateDM(da, &cda));
1762: PetscCall(DMGetCoordinatesLocal(da, &coords));
1763: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1764: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1765: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1767: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1768: PetscCall(PetscMalloc1(nx, &bc_vals));
1770: /* init the entries to -1 so VecSetValues will ignore them */
1771: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1773: j = ny - 1;
1774: for (i = 0; i < nx; i++) {
1775: PetscInt local_id;
1777: local_id = i + j * nx;
1779: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1781: bc_vals[i] = 0.0;
1782: }
1783: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1784: nbcs = 0;
1785: if ((sj + ny) == (N)) nbcs = nx;
1787: if (b) {
1788: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1789: PetscCall(VecAssemblyBegin(b));
1790: PetscCall(VecAssemblyEnd(b));
1791: }
1792: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1794: PetscCall(PetscFree(bc_vals));
1795: PetscCall(PetscFree(bc_global_ids));
1797: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1798: PetscFunctionReturn(PETSC_SUCCESS);
1799: }
1801: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
1802: {
1803: DM cda;
1804: Vec coords;
1805: PetscInt si, sj, nx, ny, i, j;
1806: PetscInt M, N;
1807: DMDACoor2d **_coords;
1808: const PetscInt *g_idx;
1809: PetscInt *bc_global_ids;
1810: PetscScalar *bc_vals;
1811: PetscInt nbcs;
1812: PetscInt n_dofs;
1813: ISLocalToGlobalMapping ltogm;
1815: PetscFunctionBeginUser;
1816: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1817: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1819: PetscCall(DMGetCoordinateDM(da, &cda));
1820: PetscCall(DMGetCoordinatesLocal(da, &coords));
1821: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1822: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1823: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1825: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1826: PetscCall(PetscMalloc1(nx, &bc_vals));
1828: /* init the entries to -1 so VecSetValues will ignore them */
1829: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1831: j = 0;
1832: for (i = 0; i < nx; i++) {
1833: PetscInt local_id;
1835: local_id = i + j * nx;
1837: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1839: bc_vals[i] = 0.0;
1840: }
1841: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1842: nbcs = 0;
1843: if (sj == 0) nbcs = nx;
1845: if (b) {
1846: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1847: PetscCall(VecAssemblyBegin(b));
1848: PetscCall(VecAssemblyEnd(b));
1849: }
1850: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1852: PetscCall(PetscFree(bc_vals));
1853: PetscCall(PetscFree(bc_global_ids));
1855: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1856: PetscFunctionReturn(PETSC_SUCCESS);
1857: }
1859: /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
1860: static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes, Mat A, Vec f)
1861: {
1862: PetscFunctionBeginUser;
1863: PetscCall(BCApplyZero_NORTH(da_Stokes, 1, A, f));
1864: PetscCall(BCApplyZero_EAST(da_Stokes, 0, A, f));
1865: PetscCall(BCApplyZero_SOUTH(da_Stokes, 1, A, f));
1866: PetscCall(BCApplyZero_WEST(da_Stokes, 0, A, f));
1867: PetscFunctionReturn(PETSC_SUCCESS);
1868: }
1870: /*TEST
1872: build:
1873: requires: !complex !single
1875: test:
1876: 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
1878: testset:
1879: 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
1880: test:
1881: suffix: 2
1882: args:
1883: output_file: output/ex43_1.out
1884: test:
1885: requires: mumps
1886: suffix: 2_mumps
1887: args: -change true -stokes_ksp_view
1888: output_file: output/ex43_2_mumps.out
1889: # mumps INFO,INFOG,RINFO,RINFOG may vary on different archs, so keep just a stable one
1890: filter: grep -E -v "(INFOG\([^7]|INFO\(|\[0\])" | grep -v " total: nonzeros="
1892: test:
1893: suffix: 3
1894: nsize: 4
1895: 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
1897: test:
1898: suffix: 4
1899: nsize: 4
1900: 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
1902: test:
1903: suffix: 5
1904: nsize: 4
1905: 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
1907: test:
1908: suffix: 6
1909: nsize: 8
1910: 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
1912: test:
1913: suffix: bjacobi
1914: nsize: 4
1915: args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type aij -stokes_ksp_converged_reason
1917: test:
1918: suffix: bjacobi_baij
1919: nsize: 4
1920: args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type baij -stokes_ksp_converged_reason
1921: output_file: output/ex43_bjacobi.out
1923: test:
1924: suffix: nested_gmg
1925: nsize: 4
1926: 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
1928: test:
1929: suffix: fetidp
1930: nsize: 8
1931: 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
1933: test:
1934: suffix: fetidp_unsym
1935: nsize: 8
1936: 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
1938: test:
1939: suffix: bddc_stokes_deluxe
1940: nsize: 8
1941: 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
1943: test:
1944: suffix: bddc_stokes_subdomainjump_deluxe
1945: nsize: 9
1946: 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
1948: TEST*/