Actual source code: ex42.c
1: static char help[] = "Solves the incompressible, variable viscosity stokes equation in 3d using Q1Q1 elements, \n\
2: stabilized with Bochev's polynomial projection method. Note that implementation here assumes \n\
3: all boundaries are free-slip, i.e. zero normal flow and zero tangential stress \n\
4: -mx : number elements in x-direction \n\
5: -my : number elements in y-direction \n\
6: -mz : number elements in z-direction \n\
7: -stokes_ksp_monitor_blocks : active monitor for each component u,v,w,p \n\
8: -model : defines viscosity and forcing function. Choose either: 0 (isoviscous), 1 (manufactured-broken), 2 (solcx-2d), 3 (sinker), 4 (subdomain jumps) \n";
10: /* Contributed by Dave May */
12: #include <petscksp.h>
13: #include <petscdmda.h>
15: #define PROFILE_TIMING
16: #define ASSEMBLE_LOWER_TRIANGULAR
18: #define NSD 3 /* number of spatial dimensions */
19: #define NODES_PER_EL 8 /* nodes per element */
20: #define U_DOFS 3 /* degrees of freedom per velocity node */
21: #define P_DOFS 1 /* degrees of freedom per pressure node */
22: #define GAUSS_POINTS 8
24: /* Gauss point based evaluation */
25: typedef struct {
26: PetscScalar gp_coords[NSD * GAUSS_POINTS];
27: PetscScalar eta[GAUSS_POINTS];
28: PetscScalar fx[GAUSS_POINTS];
29: PetscScalar fy[GAUSS_POINTS];
30: PetscScalar fz[GAUSS_POINTS];
31: PetscScalar hc[GAUSS_POINTS];
32: } GaussPointCoefficients;
34: typedef struct {
35: PetscScalar u_dof;
36: PetscScalar v_dof;
37: PetscScalar w_dof;
38: PetscScalar p_dof;
39: } StokesDOF;
41: typedef struct _p_CellProperties *CellProperties;
42: struct _p_CellProperties {
43: PetscInt ncells;
44: PetscInt mx, my, mz;
45: PetscInt sex, sey, sez;
46: GaussPointCoefficients *gpc;
47: };
49: /* elements */
50: PetscErrorCode CellPropertiesCreate(DM da_stokes, CellProperties *C)
51: {
52: CellProperties cells;
53: PetscInt mx, my, mz, sex, sey, sez;
55: PetscFunctionBeginUser;
56: PetscCall(PetscNew(&cells));
58: PetscCall(DMDAGetElementsCorners(da_stokes, &sex, &sey, &sez));
59: PetscCall(DMDAGetElementsSizes(da_stokes, &mx, &my, &mz));
60: cells->mx = mx;
61: cells->my = my;
62: cells->mz = mz;
63: cells->ncells = mx * my * mz;
64: cells->sex = sex;
65: cells->sey = sey;
66: cells->sez = sez;
68: PetscCall(PetscMalloc1(mx * my * mz, &cells->gpc));
70: *C = cells;
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
75: {
76: CellProperties cells;
78: PetscFunctionBeginUser;
79: if (!C) PetscFunctionReturn(PETSC_SUCCESS);
80: cells = *C;
81: PetscCall(PetscFree(cells->gpc));
82: PetscCall(PetscFree(cells));
83: *C = NULL;
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: PetscErrorCode CellPropertiesGetCell(CellProperties C, PetscInt II, PetscInt J, PetscInt K, GaussPointCoefficients **G)
88: {
89: PetscFunctionBeginUser;
90: *G = &C->gpc[(II - C->sex) + (J - C->sey) * C->mx + (K - C->sez) * C->mx * C->my];
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /* FEM routines */
95: /*
96: Element: Local basis function ordering
97: 1-----2
98: | |
99: | |
100: 0-----3
101: */
102: static void ShapeFunctionQ13D_Evaluate(PetscScalar _xi[], PetscScalar Ni[])
103: {
104: PetscReal xi = PetscRealPart(_xi[0]);
105: PetscReal eta = PetscRealPart(_xi[1]);
106: PetscReal zeta = PetscRealPart(_xi[2]);
108: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
109: Ni[1] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
110: Ni[2] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
111: Ni[3] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
113: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
114: Ni[5] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
115: Ni[6] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
116: Ni[7] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
117: }
119: static void ShapeFunctionQ13D_Evaluate_dxi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
120: {
121: PetscReal xi = PetscRealPart(_xi[0]);
122: PetscReal eta = PetscRealPart(_xi[1]);
123: PetscReal zeta = PetscRealPart(_xi[2]);
124: /* xi deriv */
125: GNi[0][0] = -0.125 * (1.0 - eta) * (1.0 - zeta);
126: GNi[0][1] = -0.125 * (1.0 + eta) * (1.0 - zeta);
127: GNi[0][2] = 0.125 * (1.0 + eta) * (1.0 - zeta);
128: GNi[0][3] = 0.125 * (1.0 - eta) * (1.0 - zeta);
130: GNi[0][4] = -0.125 * (1.0 - eta) * (1.0 + zeta);
131: GNi[0][5] = -0.125 * (1.0 + eta) * (1.0 + zeta);
132: GNi[0][6] = 0.125 * (1.0 + eta) * (1.0 + zeta);
133: GNi[0][7] = 0.125 * (1.0 - eta) * (1.0 + zeta);
134: /* eta deriv */
135: GNi[1][0] = -0.125 * (1.0 - xi) * (1.0 - zeta);
136: GNi[1][1] = 0.125 * (1.0 - xi) * (1.0 - zeta);
137: GNi[1][2] = 0.125 * (1.0 + xi) * (1.0 - zeta);
138: GNi[1][3] = -0.125 * (1.0 + xi) * (1.0 - zeta);
140: GNi[1][4] = -0.125 * (1.0 - xi) * (1.0 + zeta);
141: GNi[1][5] = 0.125 * (1.0 - xi) * (1.0 + zeta);
142: GNi[1][6] = 0.125 * (1.0 + xi) * (1.0 + zeta);
143: GNi[1][7] = -0.125 * (1.0 + xi) * (1.0 + zeta);
144: /* zeta deriv */
145: GNi[2][0] = -0.125 * (1.0 - xi) * (1.0 - eta);
146: GNi[2][1] = -0.125 * (1.0 - xi) * (1.0 + eta);
147: GNi[2][2] = -0.125 * (1.0 + xi) * (1.0 + eta);
148: GNi[2][3] = -0.125 * (1.0 + xi) * (1.0 - eta);
150: GNi[2][4] = 0.125 * (1.0 - xi) * (1.0 - eta);
151: GNi[2][5] = 0.125 * (1.0 - xi) * (1.0 + eta);
152: GNi[2][6] = 0.125 * (1.0 + xi) * (1.0 + eta);
153: GNi[2][7] = 0.125 * (1.0 + xi) * (1.0 - eta);
154: }
156: static void matrix_inverse_3x3(PetscScalar A[3][3], PetscScalar B[3][3])
157: {
158: PetscScalar t4, t6, t8, t10, t12, t14, t17;
160: t4 = A[2][0] * A[0][1];
161: t6 = A[2][0] * A[0][2];
162: t8 = A[1][0] * A[0][1];
163: t10 = A[1][0] * A[0][2];
164: t12 = A[0][0] * A[1][1];
165: t14 = A[0][0] * A[1][2];
166: t17 = 0.1e1 / (t4 * A[1][2] - t6 * A[1][1] - t8 * A[2][2] + t10 * A[2][1] + t12 * A[2][2] - t14 * A[2][1]);
168: B[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) * t17;
169: B[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) * t17;
170: B[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) * t17;
171: B[1][0] = -(-A[2][0] * A[1][2] + A[1][0] * A[2][2]) * t17;
172: B[1][1] = (-t6 + A[0][0] * A[2][2]) * t17;
173: B[1][2] = -(-t10 + t14) * t17;
174: B[2][0] = (-A[2][0] * A[1][1] + A[1][0] * A[2][1]) * t17;
175: B[2][1] = -(-t4 + A[0][0] * A[2][1]) * t17;
176: B[2][2] = (-t8 + t12) * t17;
177: }
179: static void ShapeFunctionQ13D_Evaluate_dx(PetscScalar GNi[][NODES_PER_EL], PetscScalar GNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
180: {
181: PetscScalar J00, J01, J02, J10, J11, J12, J20, J21, J22;
182: PetscInt n;
183: PetscScalar iJ[3][3], JJ[3][3];
185: J00 = J01 = J02 = 0.0;
186: J10 = J11 = J12 = 0.0;
187: J20 = J21 = J22 = 0.0;
188: for (n = 0; n < NODES_PER_EL; n++) {
189: PetscScalar cx = coords[NSD * n + 0];
190: PetscScalar cy = coords[NSD * n + 1];
191: PetscScalar cz = coords[NSD * n + 2];
193: /* J_ij = d(x_j) / d(xi_i) */ /* J_ij = \sum _I GNi[j][I} * x_i */
194: J00 = J00 + GNi[0][n] * cx; /* J_xx */
195: J01 = J01 + GNi[0][n] * cy; /* J_xy = dx/deta */
196: J02 = J02 + GNi[0][n] * cz; /* J_xz = dx/dzeta */
198: J10 = J10 + GNi[1][n] * cx; /* J_yx = dy/dxi */
199: J11 = J11 + GNi[1][n] * cy; /* J_yy */
200: J12 = J12 + GNi[1][n] * cz; /* J_yz */
202: J20 = J20 + GNi[2][n] * cx; /* J_zx */
203: J21 = J21 + GNi[2][n] * cy; /* J_zy */
204: J22 = J22 + GNi[2][n] * cz; /* J_zz */
205: }
207: JJ[0][0] = J00;
208: JJ[0][1] = J01;
209: JJ[0][2] = J02;
210: JJ[1][0] = J10;
211: JJ[1][1] = J11;
212: JJ[1][2] = J12;
213: JJ[2][0] = J20;
214: JJ[2][1] = J21;
215: JJ[2][2] = J22;
217: matrix_inverse_3x3(JJ, iJ);
219: *det_J = J00 * J11 * J22 - J00 * J12 * J21 - J10 * J01 * J22 + J10 * J02 * J21 + J20 * J01 * J12 - J20 * J02 * J11;
221: for (n = 0; n < NODES_PER_EL; n++) {
222: GNx[0][n] = GNi[0][n] * iJ[0][0] + GNi[1][n] * iJ[0][1] + GNi[2][n] * iJ[0][2];
223: GNx[1][n] = GNi[0][n] * iJ[1][0] + GNi[1][n] * iJ[1][1] + GNi[2][n] * iJ[1][2];
224: GNx[2][n] = GNi[0][n] * iJ[2][0] + GNi[1][n] * iJ[2][1] + GNi[2][n] * iJ[2][2];
225: }
226: }
228: static void ConstructGaussQuadrature3D(PetscInt *ngp, PetscScalar gp_xi[][NSD], PetscScalar gp_weight[])
229: {
230: *ngp = 8;
231: gp_xi[0][0] = -0.57735026919;
232: gp_xi[0][1] = -0.57735026919;
233: gp_xi[0][2] = -0.57735026919;
234: gp_xi[1][0] = -0.57735026919;
235: gp_xi[1][1] = 0.57735026919;
236: gp_xi[1][2] = -0.57735026919;
237: gp_xi[2][0] = 0.57735026919;
238: gp_xi[2][1] = 0.57735026919;
239: gp_xi[2][2] = -0.57735026919;
240: gp_xi[3][0] = 0.57735026919;
241: gp_xi[3][1] = -0.57735026919;
242: gp_xi[3][2] = -0.57735026919;
244: gp_xi[4][0] = -0.57735026919;
245: gp_xi[4][1] = -0.57735026919;
246: gp_xi[4][2] = 0.57735026919;
247: gp_xi[5][0] = -0.57735026919;
248: gp_xi[5][1] = 0.57735026919;
249: gp_xi[5][2] = 0.57735026919;
250: gp_xi[6][0] = 0.57735026919;
251: gp_xi[6][1] = 0.57735026919;
252: gp_xi[6][2] = 0.57735026919;
253: gp_xi[7][0] = 0.57735026919;
254: gp_xi[7][1] = -0.57735026919;
255: gp_xi[7][2] = 0.57735026919;
257: gp_weight[0] = 1.0;
258: gp_weight[1] = 1.0;
259: gp_weight[2] = 1.0;
260: gp_weight[3] = 1.0;
262: gp_weight[4] = 1.0;
263: gp_weight[5] = 1.0;
264: gp_weight[6] = 1.0;
265: gp_weight[7] = 1.0;
266: }
268: /*
269: i,j are the element indices
270: The unknown is a vector quantity.
271: The s[].c is used to indicate the degree of freedom.
272: */
273: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[], MatStencil s_p[], PetscInt i, PetscInt j, PetscInt k)
274: {
275: PetscInt n;
277: PetscFunctionBeginUser;
278: /* velocity */
279: n = 0;
280: /* node 0 */
281: s_u[n].i = i;
282: s_u[n].j = j;
283: s_u[n].k = k;
284: s_u[n].c = 0;
285: n++; /* Vx0 */
286: s_u[n].i = i;
287: s_u[n].j = j;
288: s_u[n].k = k;
289: s_u[n].c = 1;
290: n++; /* Vy0 */
291: s_u[n].i = i;
292: s_u[n].j = j;
293: s_u[n].k = k;
294: s_u[n].c = 2;
295: n++; /* Vz0 */
297: s_u[n].i = i;
298: s_u[n].j = j + 1;
299: s_u[n].k = k;
300: s_u[n].c = 0;
301: n++;
302: s_u[n].i = i;
303: s_u[n].j = j + 1;
304: s_u[n].k = k;
305: s_u[n].c = 1;
306: n++;
307: s_u[n].i = i;
308: s_u[n].j = j + 1;
309: s_u[n].k = k;
310: s_u[n].c = 2;
311: n++;
313: s_u[n].i = i + 1;
314: s_u[n].j = j + 1;
315: s_u[n].k = k;
316: s_u[n].c = 0;
317: n++;
318: s_u[n].i = i + 1;
319: s_u[n].j = j + 1;
320: s_u[n].k = k;
321: s_u[n].c = 1;
322: n++;
323: s_u[n].i = i + 1;
324: s_u[n].j = j + 1;
325: s_u[n].k = k;
326: s_u[n].c = 2;
327: n++;
329: s_u[n].i = i + 1;
330: s_u[n].j = j;
331: s_u[n].k = k;
332: s_u[n].c = 0;
333: n++;
334: s_u[n].i = i + 1;
335: s_u[n].j = j;
336: s_u[n].k = k;
337: s_u[n].c = 1;
338: n++;
339: s_u[n].i = i + 1;
340: s_u[n].j = j;
341: s_u[n].k = k;
342: s_u[n].c = 2;
343: n++;
345: /* */
346: s_u[n].i = i;
347: s_u[n].j = j;
348: s_u[n].k = k + 1;
349: s_u[n].c = 0;
350: n++; /* Vx4 */
351: s_u[n].i = i;
352: s_u[n].j = j;
353: s_u[n].k = k + 1;
354: s_u[n].c = 1;
355: n++; /* Vy4 */
356: s_u[n].i = i;
357: s_u[n].j = j;
358: s_u[n].k = k + 1;
359: s_u[n].c = 2;
360: n++; /* Vz4 */
362: s_u[n].i = i;
363: s_u[n].j = j + 1;
364: s_u[n].k = k + 1;
365: s_u[n].c = 0;
366: n++;
367: s_u[n].i = i;
368: s_u[n].j = j + 1;
369: s_u[n].k = k + 1;
370: s_u[n].c = 1;
371: n++;
372: s_u[n].i = i;
373: s_u[n].j = j + 1;
374: s_u[n].k = k + 1;
375: s_u[n].c = 2;
376: n++;
378: s_u[n].i = i + 1;
379: s_u[n].j = j + 1;
380: s_u[n].k = k + 1;
381: s_u[n].c = 0;
382: n++;
383: s_u[n].i = i + 1;
384: s_u[n].j = j + 1;
385: s_u[n].k = k + 1;
386: s_u[n].c = 1;
387: n++;
388: s_u[n].i = i + 1;
389: s_u[n].j = j + 1;
390: s_u[n].k = k + 1;
391: s_u[n].c = 2;
392: n++;
394: s_u[n].i = i + 1;
395: s_u[n].j = j;
396: s_u[n].k = k + 1;
397: s_u[n].c = 0;
398: n++;
399: s_u[n].i = i + 1;
400: s_u[n].j = j;
401: s_u[n].k = k + 1;
402: s_u[n].c = 1;
403: n++;
404: s_u[n].i = i + 1;
405: s_u[n].j = j;
406: s_u[n].k = k + 1;
407: s_u[n].c = 2;
408: n++;
410: /* pressure */
411: n = 0;
413: s_p[n].i = i;
414: s_p[n].j = j;
415: s_p[n].k = k;
416: s_p[n].c = 3;
417: n++; /* P0 */
418: s_p[n].i = i;
419: s_p[n].j = j + 1;
420: s_p[n].k = k;
421: s_p[n].c = 3;
422: n++;
423: s_p[n].i = i + 1;
424: s_p[n].j = j + 1;
425: s_p[n].k = k;
426: s_p[n].c = 3;
427: n++;
428: s_p[n].i = i + 1;
429: s_p[n].j = j;
430: s_p[n].k = k;
431: s_p[n].c = 3;
432: n++;
434: s_p[n].i = i;
435: s_p[n].j = j;
436: s_p[n].k = k + 1;
437: s_p[n].c = 3;
438: n++; /* P0 */
439: s_p[n].i = i;
440: s_p[n].j = j + 1;
441: s_p[n].k = k + 1;
442: s_p[n].c = 3;
443: n++;
444: s_p[n].i = i + 1;
445: s_p[n].j = j + 1;
446: s_p[n].k = k + 1;
447: s_p[n].c = 3;
448: n++;
449: s_p[n].i = i + 1;
450: s_p[n].j = j;
451: s_p[n].k = k + 1;
452: s_p[n].c = 3;
453: n++;
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords, PetscInt i, PetscInt j, PetscInt k, PetscScalar el_coord[])
458: {
459: PetscFunctionBeginUser;
460: /* get coords for the element */
461: el_coord[0] = coords[k][j][i].x;
462: el_coord[1] = coords[k][j][i].y;
463: el_coord[2] = coords[k][j][i].z;
465: el_coord[3] = coords[k][j + 1][i].x;
466: el_coord[4] = coords[k][j + 1][i].y;
467: el_coord[5] = coords[k][j + 1][i].z;
469: el_coord[6] = coords[k][j + 1][i + 1].x;
470: el_coord[7] = coords[k][j + 1][i + 1].y;
471: el_coord[8] = coords[k][j + 1][i + 1].z;
473: el_coord[9] = coords[k][j][i + 1].x;
474: el_coord[10] = coords[k][j][i + 1].y;
475: el_coord[11] = coords[k][j][i + 1].z;
477: el_coord[12] = coords[k + 1][j][i].x;
478: el_coord[13] = coords[k + 1][j][i].y;
479: el_coord[14] = coords[k + 1][j][i].z;
481: el_coord[15] = coords[k + 1][j + 1][i].x;
482: el_coord[16] = coords[k + 1][j + 1][i].y;
483: el_coord[17] = coords[k + 1][j + 1][i].z;
485: el_coord[18] = coords[k + 1][j + 1][i + 1].x;
486: el_coord[19] = coords[k + 1][j + 1][i + 1].y;
487: el_coord[20] = coords[k + 1][j + 1][i + 1].z;
489: el_coord[21] = coords[k + 1][j][i + 1].x;
490: el_coord[22] = coords[k + 1][j][i + 1].y;
491: el_coord[23] = coords[k + 1][j][i + 1].z;
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field, PetscInt i, PetscInt j, PetscInt k, StokesDOF nodal_fields[])
496: {
497: PetscFunctionBeginUser;
498: /* get the nodal fields for u */
499: nodal_fields[0].u_dof = field[k][j][i].u_dof;
500: nodal_fields[0].v_dof = field[k][j][i].v_dof;
501: nodal_fields[0].w_dof = field[k][j][i].w_dof;
503: nodal_fields[1].u_dof = field[k][j + 1][i].u_dof;
504: nodal_fields[1].v_dof = field[k][j + 1][i].v_dof;
505: nodal_fields[1].w_dof = field[k][j + 1][i].w_dof;
507: nodal_fields[2].u_dof = field[k][j + 1][i + 1].u_dof;
508: nodal_fields[2].v_dof = field[k][j + 1][i + 1].v_dof;
509: nodal_fields[2].w_dof = field[k][j + 1][i + 1].w_dof;
511: nodal_fields[3].u_dof = field[k][j][i + 1].u_dof;
512: nodal_fields[3].v_dof = field[k][j][i + 1].v_dof;
513: nodal_fields[3].w_dof = field[k][j][i + 1].w_dof;
515: nodal_fields[4].u_dof = field[k + 1][j][i].u_dof;
516: nodal_fields[4].v_dof = field[k + 1][j][i].v_dof;
517: nodal_fields[4].w_dof = field[k + 1][j][i].w_dof;
519: nodal_fields[5].u_dof = field[k + 1][j + 1][i].u_dof;
520: nodal_fields[5].v_dof = field[k + 1][j + 1][i].v_dof;
521: nodal_fields[5].w_dof = field[k + 1][j + 1][i].w_dof;
523: nodal_fields[6].u_dof = field[k + 1][j + 1][i + 1].u_dof;
524: nodal_fields[6].v_dof = field[k + 1][j + 1][i + 1].v_dof;
525: nodal_fields[6].w_dof = field[k + 1][j + 1][i + 1].w_dof;
527: nodal_fields[7].u_dof = field[k + 1][j][i + 1].u_dof;
528: nodal_fields[7].v_dof = field[k + 1][j][i + 1].v_dof;
529: nodal_fields[7].w_dof = field[k + 1][j][i + 1].w_dof;
531: /* get the nodal fields for p */
532: nodal_fields[0].p_dof = field[k][j][i].p_dof;
533: nodal_fields[1].p_dof = field[k][j + 1][i].p_dof;
534: nodal_fields[2].p_dof = field[k][j + 1][i + 1].p_dof;
535: nodal_fields[3].p_dof = field[k][j][i + 1].p_dof;
537: nodal_fields[4].p_dof = field[k + 1][j][i].p_dof;
538: nodal_fields[5].p_dof = field[k + 1][j + 1][i].p_dof;
539: nodal_fields[6].p_dof = field[k + 1][j + 1][i + 1].p_dof;
540: nodal_fields[7].p_dof = field[k + 1][j][i + 1].p_dof;
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: 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)
545: {
546: PetscInt ij;
547: PETSC_UNUSED PetscInt r, c, nr, nc;
549: nr = w_NPE * w_dof;
550: nc = u_NPE * u_dof;
552: r = w_dof * wi + wd;
553: c = u_dof * ui + ud;
555: ij = r * nc + c;
557: return ij;
558: }
560: static PetscErrorCode DMDASetValuesLocalStencil3D_ADD_VALUES(StokesDOF ***fields_F, MatStencil u_eqn[], MatStencil p_eqn[], PetscScalar Fe_u[], PetscScalar Fe_p[])
561: {
562: PetscInt n, II, J, K;
564: PetscFunctionBeginUser;
565: for (n = 0; n < NODES_PER_EL; n++) {
566: II = u_eqn[NSD * n].i;
567: J = u_eqn[NSD * n].j;
568: K = u_eqn[NSD * n].k;
570: fields_F[K][J][II].u_dof = fields_F[K][J][II].u_dof + Fe_u[NSD * n];
572: II = u_eqn[NSD * n + 1].i;
573: J = u_eqn[NSD * n + 1].j;
574: K = u_eqn[NSD * n + 1].k;
576: fields_F[K][J][II].v_dof = fields_F[K][J][II].v_dof + Fe_u[NSD * n + 1];
578: II = u_eqn[NSD * n + 2].i;
579: J = u_eqn[NSD * n + 2].j;
580: K = u_eqn[NSD * n + 2].k;
581: fields_F[K][J][II].w_dof = fields_F[K][J][II].w_dof + Fe_u[NSD * n + 2];
583: II = p_eqn[n].i;
584: J = p_eqn[n].j;
585: K = p_eqn[n].k;
587: fields_F[K][J][II].p_dof = fields_F[K][J][II].p_dof + Fe_p[n];
588: }
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: static void FormStressOperatorQ13D(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
593: {
594: PetscInt ngp;
595: PetscScalar gp_xi[GAUSS_POINTS][NSD];
596: PetscScalar gp_weight[GAUSS_POINTS];
597: PetscInt p, i, j, k;
598: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
599: PetscScalar J_p, tildeD[6];
600: PetscScalar B[6][U_DOFS * NODES_PER_EL];
601: const PetscInt nvdof = U_DOFS * NODES_PER_EL;
603: /* define quadrature rule */
604: ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);
606: /* evaluate integral */
607: for (p = 0; p < ngp; p++) {
608: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
609: ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
611: for (i = 0; i < NODES_PER_EL; i++) {
612: PetscScalar d_dx_i = GNx_p[0][i];
613: PetscScalar d_dy_i = GNx_p[1][i];
614: PetscScalar d_dz_i = GNx_p[2][i];
616: B[0][3 * i] = d_dx_i;
617: B[0][3 * i + 1] = 0.0;
618: B[0][3 * i + 2] = 0.0;
619: B[1][3 * i] = 0.0;
620: B[1][3 * i + 1] = d_dy_i;
621: B[1][3 * i + 2] = 0.0;
622: B[2][3 * i] = 0.0;
623: B[2][3 * i + 1] = 0.0;
624: B[2][3 * i + 2] = d_dz_i;
626: B[3][3 * i] = d_dy_i;
627: B[3][3 * i + 1] = d_dx_i;
628: B[3][3 * i + 2] = 0.0; /* e_xy */
629: B[4][3 * i] = d_dz_i;
630: B[4][3 * i + 1] = 0.0;
631: B[4][3 * i + 2] = d_dx_i; /* e_xz */
632: B[5][3 * i] = 0.0;
633: B[5][3 * i + 1] = d_dz_i;
634: B[5][3 * i + 2] = d_dy_i; /* e_yz */
635: }
637: tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
638: tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
639: tildeD[2] = 2.0 * gp_weight[p] * J_p * eta[p];
641: tildeD[3] = gp_weight[p] * J_p * eta[p];
642: tildeD[4] = gp_weight[p] * J_p * eta[p];
643: tildeD[5] = gp_weight[p] * J_p * eta[p];
645: /* form Bt tildeD B */
646: /*
647: Ke_ij = Bt_ik . D_kl . B_lj
648: = B_ki . D_kl . B_lj
649: = B_ki . D_kk . B_kj
650: */
651: for (i = 0; i < nvdof; i++) {
652: for (j = i; j < nvdof; j++) {
653: for (k = 0; k < 6; k++) Ke[i * nvdof + j] += B[k][i] * tildeD[k] * B[k][j];
654: }
655: }
656: }
657: /* fill lower triangular part */
658: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
659: for (i = 0; i < nvdof; i++) {
660: for (j = i; j < nvdof; j++) Ke[j * nvdof + i] = Ke[i * nvdof + j];
661: }
662: #endif
663: }
665: static void FormGradientOperatorQ13D(PetscScalar Ke[], PetscScalar coords[])
666: {
667: PetscInt ngp;
668: PetscScalar gp_xi[GAUSS_POINTS][NSD];
669: PetscScalar gp_weight[GAUSS_POINTS];
670: PetscInt p, i, j, di;
671: PetscScalar Ni_p[NODES_PER_EL];
672: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
673: PetscScalar J_p, fac;
675: /* define quadrature rule */
676: ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);
678: /* evaluate integral */
679: for (p = 0; p < ngp; p++) {
680: ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
681: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
682: ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
683: fac = gp_weight[p] * J_p;
685: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
686: for (di = 0; di < NSD; di++) { /* u dofs */
687: for (j = 0; j < NODES_PER_EL; j++) { /* p nodes, p dofs = 1 (ie no loop) */
688: PetscInt IJ;
689: IJ = ASS_MAP_wIwDI_uJuDJ(i, di, NODES_PER_EL, 3, j, 0, NODES_PER_EL, 1);
691: Ke[IJ] = Ke[IJ] - GNx_p[di][i] * Ni_p[j] * fac;
692: }
693: }
694: }
695: }
696: }
698: static void FormDivergenceOperatorQ13D(PetscScalar De[], PetscScalar coords[])
699: {
700: PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
701: PetscInt i, j;
702: PetscInt nr_g, nc_g;
704: PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));
705: FormGradientOperatorQ13D(Ge, coords);
707: nr_g = U_DOFS * NODES_PER_EL;
708: nc_g = P_DOFS * NODES_PER_EL;
710: for (i = 0; i < nr_g; i++) {
711: for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
712: }
713: }
715: static void FormStabilisationOperatorQ13D(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
716: {
717: PetscInt ngp;
718: PetscScalar gp_xi[GAUSS_POINTS][NSD];
719: PetscScalar gp_weight[GAUSS_POINTS];
720: PetscInt p, i, j;
721: PetscScalar Ni_p[NODES_PER_EL];
722: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
723: PetscScalar J_p, fac, eta_avg;
725: /* define quadrature rule */
726: ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);
728: /* evaluate integral */
729: for (p = 0; p < ngp; p++) {
730: ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
731: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
732: ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
733: fac = gp_weight[p] * J_p;
734: /*
735: for (i = 0; i < NODES_PER_EL; i++) {
736: for (j = i; j < NODES_PER_EL; j++) {
737: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
738: }
739: }
740: */
742: for (i = 0; i < NODES_PER_EL; i++) {
743: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * (Ni_p[i] * Ni_p[j] - 0.015625);
744: }
745: }
747: /* scale */
748: eta_avg = 0.0;
749: for (p = 0; p < ngp; p++) eta_avg += eta[p];
750: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
751: fac = 1.0 / eta_avg;
753: /*
754: for (i = 0; i < NODES_PER_EL; i++) {
755: for (j = i; j < NODES_PER_EL; j++) {
756: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
757: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
758: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
759: #endif
760: }
761: }
762: */
764: for (i = 0; i < NODES_PER_EL; i++) {
765: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
766: }
767: }
769: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
770: {
771: PetscInt ngp;
772: PetscScalar gp_xi[GAUSS_POINTS][NSD];
773: PetscScalar gp_weight[GAUSS_POINTS];
774: PetscInt p, i, j;
775: PetscScalar Ni_p[NODES_PER_EL];
776: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
777: PetscScalar J_p, fac, eta_avg;
779: /* define quadrature rule */
780: ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);
782: /* evaluate integral */
783: for (p = 0; p < ngp; p++) {
784: ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
785: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
786: ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
787: fac = gp_weight[p] * J_p;
789: /*
790: for (i = 0; i < NODES_PER_EL; i++) {
791: for (j = i; j < NODES_PER_EL; j++) {
792: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
793: }
794: }
795: */
797: for (i = 0; i < NODES_PER_EL; i++) {
798: 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];
799: }
800: }
802: /* scale */
803: eta_avg = 0.0;
804: for (p = 0; p < ngp; p++) eta_avg += eta[p];
805: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
806: fac = 1.0 / eta_avg;
807: /*
808: for (i = 0; i < NODES_PER_EL; i++) {
809: for (j = i; j < NODES_PER_EL; j++) {
810: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
811: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
812: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
813: #endif
814: }
815: }
816: */
818: for (i = 0; i < NODES_PER_EL; i++) {
819: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
820: }
821: }
823: static void FormMomentumRhsQ13D(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[], PetscScalar fz[])
824: {
825: PetscInt ngp;
826: PetscScalar gp_xi[GAUSS_POINTS][NSD];
827: PetscScalar gp_weight[GAUSS_POINTS];
828: PetscInt p, i;
829: PetscScalar Ni_p[NODES_PER_EL];
830: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
831: PetscScalar J_p, fac;
833: /* define quadrature rule */
834: ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);
836: /* evaluate integral */
837: for (p = 0; p < ngp; p++) {
838: ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
839: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
840: ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
841: fac = gp_weight[p] * J_p;
843: for (i = 0; i < NODES_PER_EL; i++) {
844: Fe[NSD * i] -= fac * Ni_p[i] * fx[p];
845: Fe[NSD * i + 1] -= fac * Ni_p[i] * fy[p];
846: Fe[NSD * i + 2] -= fac * Ni_p[i] * fz[p];
847: }
848: }
849: }
851: static void FormContinuityRhsQ13D(PetscScalar Fe[], PetscScalar coords[], PetscScalar hc[])
852: {
853: PetscInt ngp;
854: PetscScalar gp_xi[GAUSS_POINTS][NSD];
855: PetscScalar gp_weight[GAUSS_POINTS];
856: PetscInt p, i;
857: PetscScalar Ni_p[NODES_PER_EL];
858: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
859: PetscScalar J_p, fac;
861: /* define quadrature rule */
862: ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);
864: /* evaluate integral */
865: for (p = 0; p < ngp; p++) {
866: ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
867: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
868: ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, coords, &J_p);
869: fac = gp_weight[p] * J_p;
871: for (i = 0; i < NODES_PER_EL; i++) Fe[i] -= fac * Ni_p[i] * hc[p];
872: }
873: }
875: #define _ZERO_ROWCOL_i(A, i) \
876: { \
877: PetscInt KK; \
878: PetscScalar tmp = A[24 * (i) + (i)]; \
879: for (KK = 0; KK < 24; KK++) A[24 * (i) + KK] = 0.0; \
880: for (KK = 0; KK < 24; KK++) A[24 * KK + (i)] = 0.0; \
881: A[24 * (i) + (i)] = tmp; \
882: }
884: #define _ZERO_ROW_i(A, i) \
885: { \
886: PetscInt KK; \
887: for (KK = 0; KK < 8; KK++) A[8 * (i) + KK] = 0.0; \
888: }
890: #define _ZERO_COL_i(A, i) \
891: { \
892: PetscInt KK; \
893: for (KK = 0; KK < 8; KK++) A[24 * KK + (i)] = 0.0; \
894: }
896: static PetscErrorCode AssembleA_Stokes(Mat A, DM stokes_da, CellProperties cell_properties)
897: {
898: DM cda;
899: Vec coords;
900: DMDACoor3d ***_coords;
901: MatStencil u_eqn[NODES_PER_EL * U_DOFS];
902: MatStencil p_eqn[NODES_PER_EL * P_DOFS];
903: PetscInt sex, sey, sez, mx, my, mz;
904: PetscInt ei, ej, ek;
905: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
906: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
907: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
908: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
909: PetscScalar el_coords[NODES_PER_EL * NSD];
910: GaussPointCoefficients *props;
911: PetscScalar *prop_eta;
912: PetscInt n, M, N, P;
914: PetscFunctionBeginUser;
915: PetscCall(DMDAGetInfo(stokes_da, 0, &M, &N, &P, 0, 0, 0, 0, 0, 0, 0, 0, 0));
916: /* setup for coords */
917: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
918: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
919: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
920: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez));
921: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, &mz));
922: for (ek = sez; ek < sez + mz; ek++) {
923: for (ej = sey; ej < sey + my; ej++) {
924: for (ei = sex; ei < sex + mx; ei++) {
925: /* get coords for the element */
926: PetscCall(GetElementCoords3D(_coords, ei, ej, ek, el_coords));
927: /* get cell properties */
928: PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &props));
929: /* get coefficients for the element */
930: prop_eta = props->eta;
932: /* initialise element stiffness matrix */
933: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
934: PetscCall(PetscMemzero(Ge, sizeof(Ge)));
935: PetscCall(PetscMemzero(De, sizeof(De)));
936: PetscCall(PetscMemzero(Ce, sizeof(Ce)));
938: /* form element stiffness matrix */
939: FormStressOperatorQ13D(Ae, el_coords, prop_eta);
940: FormGradientOperatorQ13D(Ge, el_coords);
941: /*#if defined(ASSEMBLE_LOWER_TRIANGULAR)*/
942: FormDivergenceOperatorQ13D(De, el_coords);
943: /*#endif*/
944: FormStabilisationOperatorQ13D(Ce, el_coords, prop_eta);
946: /* insert element matrix into global matrix */
947: PetscCall(DMDAGetElementEqnums3D_up(u_eqn, p_eqn, ei, ej, ek));
949: for (n = 0; n < NODES_PER_EL; n++) {
950: if ((u_eqn[3 * n].i == 0) || (u_eqn[3 * n].i == M - 1)) {
951: _ZERO_ROWCOL_i(Ae, 3 * n);
952: _ZERO_ROW_i(Ge, 3 * n);
953: _ZERO_COL_i(De, 3 * n);
954: }
956: if ((u_eqn[3 * n + 1].j == 0) || (u_eqn[3 * n + 1].j == N - 1)) {
957: _ZERO_ROWCOL_i(Ae, 3 * n + 1);
958: _ZERO_ROW_i(Ge, 3 * n + 1);
959: _ZERO_COL_i(De, 3 * n + 1);
960: }
962: if ((u_eqn[3 * n + 2].k == 0) || (u_eqn[3 * n + 2].k == P - 1)) {
963: _ZERO_ROWCOL_i(Ae, 3 * n + 2);
964: _ZERO_ROW_i(Ge, 3 * n + 2);
965: _ZERO_COL_i(De, 3 * n + 2);
966: }
967: }
968: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
969: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
970: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
971: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
972: }
973: }
974: }
975: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
976: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
978: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
980: PetscFunctionReturn(PETSC_SUCCESS);
981: }
983: static PetscErrorCode AssembleA_PCStokes(Mat A, DM stokes_da, CellProperties cell_properties)
984: {
985: DM cda;
986: Vec coords;
987: DMDACoor3d ***_coords;
988: MatStencil u_eqn[NODES_PER_EL * U_DOFS];
989: MatStencil p_eqn[NODES_PER_EL * P_DOFS];
990: PetscInt sex, sey, sez, mx, my, mz;
991: PetscInt ei, ej, ek;
992: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
993: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
994: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
995: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
996: PetscScalar el_coords[NODES_PER_EL * NSD];
997: GaussPointCoefficients *props;
998: PetscScalar *prop_eta;
999: PetscInt n, M, N, P;
1001: PetscFunctionBeginUser;
1002: PetscCall(DMDAGetInfo(stokes_da, 0, &M, &N, &P, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1003: /* setup for coords */
1004: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
1005: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
1006: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1008: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez));
1009: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, &mz));
1010: for (ek = sez; ek < sez + mz; ek++) {
1011: for (ej = sey; ej < sey + my; ej++) {
1012: for (ei = sex; ei < sex + mx; ei++) {
1013: /* get coords for the element */
1014: PetscCall(GetElementCoords3D(_coords, ei, ej, ek, el_coords));
1015: /* get cell properties */
1016: PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &props));
1017: /* get coefficients for the element */
1018: prop_eta = props->eta;
1020: /* initialise element stiffness matrix */
1021: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
1022: PetscCall(PetscMemzero(Ge, sizeof(Ge)));
1023: PetscCall(PetscMemzero(De, sizeof(De)));
1024: PetscCall(PetscMemzero(Ce, sizeof(Ce)));
1026: /* form element stiffness matrix */
1027: FormStressOperatorQ13D(Ae, el_coords, prop_eta);
1028: FormGradientOperatorQ13D(Ge, el_coords);
1029: /* FormDivergenceOperatorQ13D(De,el_coords); */
1030: FormScaledMassMatrixOperatorQ13D(Ce, el_coords, prop_eta);
1032: /* insert element matrix into global matrix */
1033: PetscCall(DMDAGetElementEqnums3D_up(u_eqn, p_eqn, ei, ej, ek));
1035: for (n = 0; n < NODES_PER_EL; n++) {
1036: if ((u_eqn[3 * n].i == 0) || (u_eqn[3 * n].i == M - 1)) {
1037: _ZERO_ROWCOL_i(Ae, 3 * n);
1038: _ZERO_ROW_i(Ge, 3 * n);
1039: _ZERO_COL_i(De, 3 * n);
1040: }
1042: if ((u_eqn[3 * n + 1].j == 0) || (u_eqn[3 * n + 1].j == N - 1)) {
1043: _ZERO_ROWCOL_i(Ae, 3 * n + 1);
1044: _ZERO_ROW_i(Ge, 3 * n + 1);
1045: _ZERO_COL_i(De, 3 * n + 1);
1046: }
1048: if ((u_eqn[3 * n + 2].k == 0) || (u_eqn[3 * n + 2].k == P - 1)) {
1049: _ZERO_ROWCOL_i(Ae, 3 * n + 2);
1050: _ZERO_ROW_i(Ge, 3 * n + 2);
1051: _ZERO_COL_i(De, 3 * n + 2);
1052: }
1053: }
1054: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
1055: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
1056: /*PetscCall(MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES));*/
1057: PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
1058: }
1059: }
1060: }
1061: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1062: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1064: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1065: PetscFunctionReturn(PETSC_SUCCESS);
1066: }
1068: static PetscErrorCode AssembleF_Stokes(Vec F, DM stokes_da, CellProperties cell_properties)
1069: {
1070: DM cda;
1071: Vec coords;
1072: DMDACoor3d ***_coords;
1073: MatStencil u_eqn[NODES_PER_EL * U_DOFS];
1074: MatStencil p_eqn[NODES_PER_EL * P_DOFS];
1075: PetscInt sex, sey, sez, mx, my, mz;
1076: PetscInt ei, ej, ek;
1077: PetscScalar Fe[NODES_PER_EL * U_DOFS];
1078: PetscScalar He[NODES_PER_EL * P_DOFS];
1079: PetscScalar el_coords[NODES_PER_EL * NSD];
1080: GaussPointCoefficients *props;
1081: PetscScalar *prop_fx, *prop_fy, *prop_fz, *prop_hc;
1082: Vec local_F;
1083: StokesDOF ***ff;
1084: PetscInt n, M, N, P;
1086: PetscFunctionBeginUser;
1087: PetscCall(DMDAGetInfo(stokes_da, 0, &M, &N, &P, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1088: /* setup for coords */
1089: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
1090: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
1091: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1093: /* get access to the vector */
1094: PetscCall(DMGetLocalVector(stokes_da, &local_F));
1095: PetscCall(VecZeroEntries(local_F));
1096: PetscCall(DMDAVecGetArray(stokes_da, local_F, &ff));
1097: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez));
1098: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, &mz));
1099: for (ek = sez; ek < sez + mz; ek++) {
1100: for (ej = sey; ej < sey + my; ej++) {
1101: for (ei = sex; ei < sex + mx; ei++) {
1102: /* get coords for the element */
1103: PetscCall(GetElementCoords3D(_coords, ei, ej, ek, el_coords));
1104: /* get cell properties */
1105: PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &props));
1106: /* get coefficients for the element */
1107: prop_fx = props->fx;
1108: prop_fy = props->fy;
1109: prop_fz = props->fz;
1110: prop_hc = props->hc;
1112: /* initialise element stiffness matrix */
1113: PetscCall(PetscMemzero(Fe, sizeof(Fe)));
1114: PetscCall(PetscMemzero(He, sizeof(He)));
1116: /* form element stiffness matrix */
1117: FormMomentumRhsQ13D(Fe, el_coords, prop_fx, prop_fy, prop_fz);
1118: FormContinuityRhsQ13D(He, el_coords, prop_hc);
1120: /* insert element matrix into global matrix */
1121: PetscCall(DMDAGetElementEqnums3D_up(u_eqn, p_eqn, ei, ej, ek));
1123: for (n = 0; n < NODES_PER_EL; n++) {
1124: if ((u_eqn[3 * n].i == 0) || (u_eqn[3 * n].i == M - 1)) Fe[3 * n] = 0.0;
1126: if ((u_eqn[3 * n + 1].j == 0) || (u_eqn[3 * n + 1].j == N - 1)) Fe[3 * n + 1] = 0.0;
1128: if ((u_eqn[3 * n + 2].k == 0) || (u_eqn[3 * n + 2].k == P - 1)) Fe[3 * n + 2] = 0.0;
1129: }
1131: PetscCall(DMDASetValuesLocalStencil3D_ADD_VALUES(ff, u_eqn, p_eqn, Fe, He));
1132: }
1133: }
1134: }
1135: PetscCall(DMDAVecRestoreArray(stokes_da, local_F, &ff));
1136: PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
1137: PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
1138: PetscCall(DMRestoreLocalVector(stokes_da, &local_F));
1140: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1141: PetscFunctionReturn(PETSC_SUCCESS);
1142: }
1144: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta, PetscReal *MX, PetscReal *MY, PetscReal *MZ)
1145: {
1146: *theta = 0.0;
1147: *MX = 2.0 * PETSC_PI;
1148: *MY = 2.0 * PETSC_PI;
1149: *MZ = 2.0 * PETSC_PI;
1150: }
1151: static void evaluate_MS_FrankKamentski(PetscReal pos[], PetscReal v[], PetscReal *p, PetscReal *eta, PetscReal Fm[], PetscReal *Fc)
1152: {
1153: PetscReal x, y, z;
1154: PetscReal theta, MX, MY, MZ;
1156: evaluate_MS_FrankKamentski_constants(&theta, &MX, &MY, &MZ);
1157: x = pos[0];
1158: y = pos[1];
1159: z = pos[2];
1160: if (v) {
1161: /*
1162: v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1163: v[1] = z*cos(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1164: v[2] = -(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z);
1165: */
1166: /*
1167: v[0] = PetscSinReal(PETSC_PI*x);
1168: v[1] = PetscSinReal(PETSC_PI*y);
1169: v[2] = PetscSinReal(PETSC_PI*z);
1170: */
1171: v[0] = PetscPowRealInt(z, 3) * PetscExpReal(y) * PetscSinReal(PETSC_PI * x);
1172: v[1] = z * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y);
1173: v[2] = PetscPowRealInt(z, 2) * (PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 - PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) / 2) - PETSC_PI * PetscPowRealInt(z, 4) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) / 4;
1174: }
1175: if (p) *p = PetscPowRealInt(x, 2) + PetscPowRealInt(y, 2) + PetscPowRealInt(z, 2);
1176: if (eta) {
1177: /* eta = PetscExpReal(-theta*(1.0 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));*/
1178: *eta = 1.0;
1179: }
1180: if (Fm) {
1181: /*
1182: Fm[0] = -2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(PETSC_PI*x) - 0.2*MZ*theta*(-1.5*PetscPowRealInt(x,2)*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PetscPowRealInt(z,2)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x))*PetscCosReal(MX*x)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) - 0.2*PETSC_PI*MX*theta*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscCosReal(MZ*z)*PetscExpReal(y)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) + 2.0*(3.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 3.0*PETSC_PI*PetscPowRealInt(x,2)*PetscCosReal(2.0*PETSC_PI*z))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 1.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) ;
1183: Fm[1] = -2*y - 0.2*MX*theta*(0.5*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 1.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x))*PetscCosReal(MZ*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) - 0.2*MZ*theta*(-1.5*PetscPowRealInt(y,2)*PetscSinReal(2.0*PETSC_PI*z) + 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y))*PetscCosReal(MX*x)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) + 2.0*(-2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.5*PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 2*PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(-z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) - 6.0*PETSC_PI*PetscPowRealInt(y,2)*PetscCosReal(2.0*PETSC_PI*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));
1184: Fm[2] = -2*z + 8.0*PetscPowRealInt(PETSC_PI,2)*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(2.0*PETSC_PI*z) - 0.2*MX*theta*(-1.5*PetscPowRealInt(x,2)*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PetscPowRealInt(z,2)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x))*PetscCosReal(MZ*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MX*x)*PetscSinReal(MY*y) + 0.4*PETSC_PI*MZ*theta*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(MX*x)*PetscCosReal(2.0*PETSC_PI*z)*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)))*PetscSinReal(MY*y)*PetscSinReal(MZ*z) + 2.0*(-3.0*x*PetscSinReal(2.0*PETSC_PI*z) + 1.5*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*(-3.0*y*PetscSinReal(2.0*PETSC_PI*z) - 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 0.5*PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) + 2.0*theta*(1 + 0.1*MY*PetscCosReal(MX*x)*PetscCosReal(MY*y)*PetscCosReal(MZ*z))*(-1.5*PetscPowRealInt(y,2)*PetscSinReal(2.0*PETSC_PI*z) + 0.5*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y))*PetscExpReal(-theta*(1 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y))) ;
1185: */
1186: /*
1187: Fm[0]=-2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*x);
1188: Fm[1]=-2*y - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*y);
1189: Fm[2]=-2*z - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*z);
1190: */
1191: /*
1192: Fm[0] = -2*x + PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 6.0*z*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) - 6.0*PETSC_PI*PetscPowRealInt(x,2)*PetscCosReal(2.0*PETSC_PI*z) - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x) + 2.0*PETSC_PI*z*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y)*PetscSinReal(2.0*PETSC_PI*x) - 2.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(PETSC_PI*y)*PetscExpReal(-y)*PetscSinReal(2.0*PETSC_PI*x) ;
1193: Fm[1] = -2*y - 6.0*PETSC_PI*PetscPowRealInt(y,2)*PetscCosReal(2.0*PETSC_PI*z) + 2.0*z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) - 6.0*z*PetscPowRealInt(PETSC_PI,2)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) - 4.0*PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);
1194: Fm[2] = -2*z - 6.0*x*PetscSinReal(2.0*PETSC_PI*z) - 6.0*y*PetscSinReal(2.0*PETSC_PI*z) - PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + 8.0*PetscPowRealInt(PETSC_PI,2)*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z) + 3.0*PETSC_PI*PetscPowRealInt(z,2)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y) ;
1195: */
1196: Fm[0] = -2 * x + 2 * z * (PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(PETSC_PI * y) * PetscExpReal(-y) * PetscSinReal(2.0 * PETSC_PI * x) - 1.0 * PETSC_PI * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) * PetscSinReal(2.0 * PETSC_PI * x)) + PetscPowRealInt(z, 3) * PetscExpReal(y) * PetscSinReal(PETSC_PI * x) + 6.0 * z * PetscExpReal(y) * PetscSinReal(PETSC_PI * x) - 1.0 * PetscPowRealInt(PETSC_PI, 2) * PetscPowRealInt(z, 3) * PetscExpReal(y) * PetscSinReal(PETSC_PI * x) + 2.0 * PETSC_PI * z * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) * PetscSinReal(2.0 * PETSC_PI * x) - 2.0 * z * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(PETSC_PI * y) * PetscExpReal(-y) * PetscSinReal(2.0 * PETSC_PI * x);
1197: Fm[1] = -2 * y + 2 * z * (-PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 + PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 + PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y)) + 2.0 * z * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) - 6.0 * z * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) - 4.0 * PETSC_PI * z * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y);
1198: Fm[2] = -2 * z + PetscPowRealInt(z, 2) * (-2.0 * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) + 2.0 * PetscPowRealInt(PETSC_PI, 3) * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y)) + PetscPowRealInt(z, 2) * (PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 - 3 * PetscPowRealInt(PETSC_PI, 2) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) / 2 + PetscPowRealInt(PETSC_PI, 3) * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) / 2 - 3 * PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) / 2) + 1.0 * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y) * PetscSinReal(PETSC_PI * y) + 0.25 * PetscPowRealInt(PETSC_PI, 3) * PetscPowRealInt(z, 4) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) - 0.25 * PETSC_PI * PetscPowRealInt(z, 4) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) - 3.0 * PETSC_PI * PetscPowRealInt(z, 2) * PetscCosReal(PETSC_PI * x) * PetscExpReal(y) - 1.0 * PETSC_PI * PetscCosReal(PETSC_PI * y) * PetscCosReal(2.0 * PETSC_PI * x) * PetscExpReal(-y);
1199: }
1200: if (Fc) {
1201: /* Fc = -2.0*PETSC_PI*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(2.0*PETSC_PI*z) - z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y) ;*/
1202: /* Fc = PETSC_PI*PetscCosReal(PETSC_PI*x) + PETSC_PI*PetscCosReal(PETSC_PI*y) + PETSC_PI*PetscCosReal(PETSC_PI*z);*/
1203: /* Fc = -2.0*PETSC_PI*(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscCosReal(2.0*PETSC_PI*z) - z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y) + PETSC_PI*PetscPowRealInt(z,3)*PetscCosReal(PETSC_PI*x)*PetscExpReal(y) + PETSC_PI*z*PetscCosReal(PETSC_PI*y)*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y);*/
1204: *Fc = 0.0;
1205: }
1206: }
1208: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx, PetscInt my, PetscInt mz, DM *_da, Vec *_X)
1209: {
1210: DM da, cda;
1211: Vec X;
1212: StokesDOF ***_stokes;
1213: Vec coords;
1214: DMDACoor3d ***_coords;
1215: PetscInt si, sj, sk, ei, ej, ek, i, j, k;
1217: PetscFunctionBeginUser;
1218: PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, mz + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 4, 1, NULL, NULL, NULL, &da));
1219: PetscCall(DMSetFromOptions(da));
1220: PetscCall(DMSetUp(da));
1221: PetscCall(DMDASetFieldName(da, 0, "anlytic_Vx"));
1222: PetscCall(DMDASetFieldName(da, 1, "anlytic_Vy"));
1223: PetscCall(DMDASetFieldName(da, 2, "anlytic_Vz"));
1224: PetscCall(DMDASetFieldName(da, 3, "analytic_P"));
1226: PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
1228: PetscCall(DMGetCoordinatesLocal(da, &coords));
1229: PetscCall(DMGetCoordinateDM(da, &cda));
1230: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1232: PetscCall(DMCreateGlobalVector(da, &X));
1233: PetscCall(DMDAVecGetArray(da, X, &_stokes));
1235: PetscCall(DMDAGetCorners(da, &si, &sj, &sk, &ei, &ej, &ek));
1236: for (k = sk; k < sk + ek; k++) {
1237: for (j = sj; j < sj + ej; j++) {
1238: for (i = si; i < si + ei; i++) {
1239: PetscReal pos[NSD], pressure, vel[NSD];
1241: pos[0] = PetscRealPart(_coords[k][j][i].x);
1242: pos[1] = PetscRealPart(_coords[k][j][i].y);
1243: pos[2] = PetscRealPart(_coords[k][j][i].z);
1245: evaluate_MS_FrankKamentski(pos, vel, &pressure, NULL, NULL, NULL);
1247: _stokes[k][j][i].u_dof = vel[0];
1248: _stokes[k][j][i].v_dof = vel[1];
1249: _stokes[k][j][i].w_dof = vel[2];
1250: _stokes[k][j][i].p_dof = pressure;
1251: }
1252: }
1253: }
1254: PetscCall(DMDAVecRestoreArray(da, X, &_stokes));
1255: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1257: *_da = da;
1258: *_X = X;
1259: PetscFunctionReturn(PETSC_SUCCESS);
1260: }
1262: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da, Vec X, Vec X_analytic)
1263: {
1264: DM cda;
1265: Vec coords, X_analytic_local, X_local;
1266: DMDACoor3d ***_coords;
1267: PetscInt sex, sey, sez, mx, my, mz;
1268: PetscInt ei, ej, ek;
1269: PetscScalar el_coords[NODES_PER_EL * NSD];
1270: StokesDOF ***stokes_analytic, ***stokes;
1271: StokesDOF stokes_analytic_e[NODES_PER_EL], stokes_e[NODES_PER_EL];
1273: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
1274: PetscScalar Ni_p[NODES_PER_EL];
1275: PetscInt ngp;
1276: PetscScalar gp_xi[GAUSS_POINTS][NSD];
1277: PetscScalar gp_weight[GAUSS_POINTS];
1278: PetscInt p, i;
1279: PetscScalar J_p, fac;
1280: PetscScalar h, p_e_L2, u_e_L2, u_e_H1, p_L2, u_L2, u_H1, tp_L2, tu_L2, tu_H1;
1281: PetscScalar tint_p_ms, tint_p, int_p_ms, int_p;
1282: PetscInt M;
1283: PetscReal xymin[NSD], xymax[NSD];
1285: PetscFunctionBeginUser;
1286: /* define quadrature rule */
1287: ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);
1289: /* setup for coords */
1290: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
1291: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
1292: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1294: /* setup for analytic */
1295: PetscCall(DMCreateLocalVector(stokes_da, &X_analytic_local));
1296: PetscCall(DMGlobalToLocalBegin(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
1297: PetscCall(DMGlobalToLocalEnd(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
1298: PetscCall(DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic));
1300: /* setup for solution */
1301: PetscCall(DMCreateLocalVector(stokes_da, &X_local));
1302: PetscCall(DMGlobalToLocalBegin(stokes_da, X, INSERT_VALUES, X_local));
1303: PetscCall(DMGlobalToLocalEnd(stokes_da, X, INSERT_VALUES, X_local));
1304: PetscCall(DMDAVecGetArray(stokes_da, X_local, &stokes));
1306: PetscCall(DMDAGetInfo(stokes_da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1307: PetscCall(DMGetBoundingBox(stokes_da, xymin, xymax));
1309: h = (xymax[0] - xymin[0]) / ((PetscReal)(M - 1));
1311: tp_L2 = tu_L2 = tu_H1 = 0.0;
1312: tint_p_ms = tint_p = 0.0;
1314: PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, &sez));
1315: PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, &mz));
1316: for (ek = sez; ek < sez + mz; ek++) {
1317: for (ej = sey; ej < sey + my; ej++) {
1318: for (ei = sex; ei < sex + mx; ei++) {
1319: /* get coords for the element */
1320: PetscCall(GetElementCoords3D(_coords, ei, ej, ek, el_coords));
1321: PetscCall(StokesDAGetNodalFields3D(stokes, ei, ej, ek, stokes_e));
1322: PetscCall(StokesDAGetNodalFields3D(stokes_analytic, ei, ej, ek, stokes_analytic_e));
1324: /* evaluate integral */
1325: for (p = 0; p < ngp; p++) {
1326: ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
1327: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
1328: ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, el_coords, &J_p);
1329: fac = gp_weight[p] * J_p;
1331: for (i = 0; i < NODES_PER_EL; i++) {
1332: tint_p_ms = tint_p_ms + fac * Ni_p[i] * stokes_analytic_e[i].p_dof;
1333: tint_p = tint_p + fac * Ni_p[i] * stokes_e[i].p_dof;
1334: }
1335: }
1336: }
1337: }
1338: }
1339: PetscCall(MPIU_Allreduce(&tint_p_ms, &int_p_ms, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1340: PetscCall(MPIU_Allreduce(&tint_p, &int_p, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1342: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\\int P dv %1.4e (h) %1.4e (ms)\n", (double)PetscRealPart(int_p), (double)PetscRealPart(int_p_ms)));
1344: /* remove mine and add manufacture one */
1345: PetscCall(DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic));
1346: PetscCall(DMDAVecRestoreArray(stokes_da, X_local, &stokes));
1348: {
1349: PetscInt k, L, dof;
1350: PetscScalar *fields;
1351: PetscCall(DMDAGetInfo(stokes_da, 0, 0, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
1353: PetscCall(VecGetLocalSize(X_local, &L));
1354: PetscCall(VecGetArray(X_local, &fields));
1355: for (k = 0; k < L / dof; k++) fields[dof * k + 3] = fields[dof * k + 3] - int_p + int_p_ms;
1356: PetscCall(VecRestoreArray(X_local, &fields));
1358: PetscCall(VecGetLocalSize(X, &L));
1359: PetscCall(VecGetArray(X, &fields));
1360: for (k = 0; k < L / dof; k++) fields[dof * k + 3] = fields[dof * k + 3] - int_p + int_p_ms;
1361: PetscCall(VecRestoreArray(X, &fields));
1362: }
1364: PetscCall(DMDAVecGetArray(stokes_da, X_local, &stokes));
1365: PetscCall(DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic));
1367: for (ek = sez; ek < sez + mz; ek++) {
1368: for (ej = sey; ej < sey + my; ej++) {
1369: for (ei = sex; ei < sex + mx; ei++) {
1370: /* get coords for the element */
1371: PetscCall(GetElementCoords3D(_coords, ei, ej, ek, el_coords));
1372: PetscCall(StokesDAGetNodalFields3D(stokes, ei, ej, ek, stokes_e));
1373: PetscCall(StokesDAGetNodalFields3D(stokes_analytic, ei, ej, ek, stokes_analytic_e));
1375: /* evaluate integral */
1376: p_e_L2 = 0.0;
1377: u_e_L2 = 0.0;
1378: u_e_H1 = 0.0;
1379: for (p = 0; p < ngp; p++) {
1380: ShapeFunctionQ13D_Evaluate(gp_xi[p], Ni_p);
1381: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p], GNi_p);
1382: ShapeFunctionQ13D_Evaluate_dx(GNi_p, GNx_p, el_coords, &J_p);
1383: fac = gp_weight[p] * J_p;
1385: for (i = 0; i < NODES_PER_EL; i++) {
1386: PetscScalar u_error, v_error, w_error;
1388: 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);
1390: u_error = stokes_e[i].u_dof - stokes_analytic_e[i].u_dof;
1391: v_error = stokes_e[i].v_dof - stokes_analytic_e[i].v_dof;
1392: w_error = stokes_e[i].w_dof - stokes_analytic_e[i].w_dof;
1393: u_e_L2 += fac * Ni_p[i] * (u_error * u_error + v_error * v_error + w_error * w_error);
1395: u_e_H1 = u_e_H1 + fac * (GNx_p[0][i] * u_error * GNx_p[0][i] * u_error /* du/dx */
1396: + GNx_p[1][i] * u_error * GNx_p[1][i] * u_error + GNx_p[2][i] * u_error * GNx_p[2][i] * u_error + GNx_p[0][i] * v_error * GNx_p[0][i] * v_error /* dv/dx */
1397: + GNx_p[1][i] * v_error * GNx_p[1][i] * v_error + GNx_p[2][i] * v_error * GNx_p[2][i] * v_error + GNx_p[0][i] * w_error * GNx_p[0][i] * w_error /* dw/dx */
1398: + GNx_p[1][i] * w_error * GNx_p[1][i] * w_error + GNx_p[2][i] * w_error * GNx_p[2][i] * w_error);
1399: }
1400: }
1402: tp_L2 += p_e_L2;
1403: tu_L2 += u_e_L2;
1404: tu_H1 += u_e_H1;
1405: }
1406: }
1407: }
1408: PetscCall(MPIU_Allreduce(&tp_L2, &p_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1409: PetscCall(MPIU_Allreduce(&tu_L2, &u_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1410: PetscCall(MPIU_Allreduce(&tu_H1, &u_H1, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1411: p_L2 = PetscSqrtScalar(p_L2);
1412: u_L2 = PetscSqrtScalar(u_L2);
1413: u_H1 = PetscSqrtScalar(u_H1);
1415: 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)));
1417: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1419: PetscCall(DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic));
1420: PetscCall(VecDestroy(&X_analytic_local));
1421: PetscCall(DMDAVecRestoreArray(stokes_da, X_local, &stokes));
1422: PetscCall(VecDestroy(&X_local));
1423: PetscFunctionReturn(PETSC_SUCCESS);
1424: }
1426: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da, Vec FIELD, const char file_prefix[])
1427: {
1428: char vtk_filename[PETSC_MAX_PATH_LEN];
1429: PetscMPIInt rank;
1430: MPI_Comm comm;
1431: FILE *vtk_fp = NULL;
1432: const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
1433: PetscInt si, sj, sk, nx, ny, nz, i;
1434: PetscInt f, n_fields, N;
1435: DM cda;
1436: Vec coords;
1437: Vec l_FIELD;
1438: PetscScalar *_L_FIELD;
1439: PetscInt memory_offset;
1440: PetscScalar *buffer;
1442: PetscFunctionBeginUser;
1444: /* create file name */
1445: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
1446: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1447: PetscCall(PetscSNPrintf(vtk_filename, sizeof(vtk_filename), "subdomain-%s-p%1.4d.vts", file_prefix, rank));
1449: /* open file and write header */
1450: vtk_fp = fopen(vtk_filename, "w");
1451: PetscCheck(vtk_fp, PETSC_COMM_SELF, PETSC_ERR_SYS, "Cannot open file = %s ", vtk_filename);
1453: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<?xml version=\"1.0\"?>\n"));
1455: /* coords */
1456: PetscCall(DMDAGetGhostCorners(da, &si, &sj, &sk, &nx, &ny, &nz));
1457: N = nx * ny * nz;
1459: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order));
1460: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <StructuredGrid WholeExtent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", si, si + nx - 1, sj, sj + ny - 1, sk, sk + nz - 1));
1461: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\">\n", si, si + nx - 1, sj, sj + ny - 1, sk, sk + nz - 1));
1463: memory_offset = 0;
1465: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <CellData></CellData>\n"));
1467: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <Points>\n"));
1469: /* copy coordinates */
1470: PetscCall(DMGetCoordinateDM(da, &cda));
1471: PetscCall(DMGetCoordinatesLocal(da, &coords));
1472: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n", memory_offset));
1473: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar) * N * 3;
1475: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </Points>\n"));
1477: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PointData Scalars=\" "));
1478: PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_fields, 0, 0, 0, 0, 0));
1479: for (f = 0; f < n_fields; f++) {
1480: const char *field_name;
1481: PetscCall(DMDAGetFieldName(da, f, &field_name));
1482: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "%s ", field_name));
1483: }
1484: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "\">\n"));
1486: for (f = 0; f < n_fields; f++) {
1487: const char *field_name;
1489: PetscCall(DMDAGetFieldName(da, f, &field_name));
1490: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <DataArray type=\"Float64\" Name=\"%s\" format=\"appended\" offset=\"%" PetscInt_FMT "\"/>\n", field_name, memory_offset));
1491: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar) * N;
1492: }
1494: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </PointData>\n"));
1495: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </Piece>\n"));
1496: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </StructuredGrid>\n"));
1498: PetscCall(PetscMalloc1(N, &buffer));
1499: PetscCall(DMGetLocalVector(da, &l_FIELD));
1500: PetscCall(DMGlobalToLocalBegin(da, FIELD, INSERT_VALUES, l_FIELD));
1501: PetscCall(DMGlobalToLocalEnd(da, FIELD, INSERT_VALUES, l_FIELD));
1502: PetscCall(VecGetArray(l_FIELD, &_L_FIELD));
1504: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <AppendedData encoding=\"raw\">\n"));
1505: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "_"));
1507: /* write coordinates */
1508: {
1509: int length = sizeof(PetscScalar) * N * 3;
1510: PetscScalar *allcoords;
1512: fwrite(&length, sizeof(int), 1, vtk_fp);
1513: PetscCall(VecGetArray(coords, &allcoords));
1514: fwrite(allcoords, sizeof(PetscScalar), 3 * N, vtk_fp);
1515: PetscCall(VecRestoreArray(coords, &allcoords));
1516: }
1517: /* write fields */
1518: for (f = 0; f < n_fields; f++) {
1519: int length = sizeof(PetscScalar) * N;
1520: fwrite(&length, sizeof(int), 1, vtk_fp);
1521: /* load */
1522: for (i = 0; i < N; i++) buffer[i] = _L_FIELD[n_fields * i + f];
1524: /* write */
1525: fwrite(buffer, sizeof(PetscScalar), N, vtk_fp);
1526: }
1527: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "\n </AppendedData>\n"));
1529: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "</VTKFile>\n"));
1531: PetscCall(PetscFree(buffer));
1532: PetscCall(VecRestoreArray(l_FIELD, &_L_FIELD));
1533: PetscCall(DMRestoreLocalVector(da, &l_FIELD));
1535: if (vtk_fp) {
1536: fclose(vtk_fp);
1537: vtk_fp = NULL;
1538: }
1540: PetscFunctionReturn(PETSC_SUCCESS);
1541: }
1543: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp, PetscInt indent_level, DM da, const char local_file_prefix[])
1544: {
1545: PetscMPIInt size, rank;
1546: MPI_Comm comm;
1547: const PetscInt *lx, *ly, *lz;
1548: PetscInt M, N, P, pM, pN, pP, sum, *olx, *oly, *olz;
1549: PetscInt *osx, *osy, *osz, *oex, *oey, *oez;
1550: PetscInt i, j, k, II, stencil;
1552: PetscFunctionBeginUser;
1553: /* create file name */
1554: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
1555: PetscCallMPI(MPI_Comm_size(comm, &size));
1556: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1558: PetscCall(DMDAGetInfo(da, 0, &M, &N, &P, &pM, &pN, &pP, 0, &stencil, 0, 0, 0, 0));
1559: PetscCall(DMDAGetOwnershipRanges(da, &lx, &ly, &lz));
1561: /* generate start,end list */
1562: PetscCall(PetscMalloc1(pM + 1, &olx));
1563: PetscCall(PetscMalloc1(pN + 1, &oly));
1564: PetscCall(PetscMalloc1(pP + 1, &olz));
1565: sum = 0;
1566: for (i = 0; i < pM; i++) {
1567: olx[i] = sum;
1568: sum = sum + lx[i];
1569: }
1570: olx[pM] = sum;
1571: sum = 0;
1572: for (i = 0; i < pN; i++) {
1573: oly[i] = sum;
1574: sum = sum + ly[i];
1575: }
1576: oly[pN] = sum;
1577: sum = 0;
1578: for (i = 0; i < pP; i++) {
1579: olz[i] = sum;
1580: sum = sum + lz[i];
1581: }
1582: olz[pP] = sum;
1584: PetscCall(PetscMalloc1(pM, &osx));
1585: PetscCall(PetscMalloc1(pN, &osy));
1586: PetscCall(PetscMalloc1(pP, &osz));
1587: PetscCall(PetscMalloc1(pM, &oex));
1588: PetscCall(PetscMalloc1(pN, &oey));
1589: PetscCall(PetscMalloc1(pP, &oez));
1590: for (i = 0; i < pM; i++) {
1591: osx[i] = olx[i] - stencil;
1592: oex[i] = olx[i] + lx[i] + stencil;
1593: if (osx[i] < 0) osx[i] = 0;
1594: if (oex[i] > M) oex[i] = M;
1595: }
1597: for (i = 0; i < pN; i++) {
1598: osy[i] = oly[i] - stencil;
1599: oey[i] = oly[i] + ly[i] + stencil;
1600: if (osy[i] < 0) osy[i] = 0;
1601: if (oey[i] > M) oey[i] = N;
1602: }
1603: for (i = 0; i < pP; i++) {
1604: osz[i] = olz[i] - stencil;
1605: oez[i] = olz[i] + lz[i] + stencil;
1606: if (osz[i] < 0) osz[i] = 0;
1607: if (oez[i] > P) oez[i] = P;
1608: }
1610: for (k = 0; k < pP; k++) {
1611: for (j = 0; j < pN; j++) {
1612: for (i = 0; i < pM; i++) {
1613: char name[PETSC_MAX_PATH_LEN];
1614: PetscInt procid = i + j * pM + k * pM * pN; /* convert proc(i,j,k) to pid */
1615: PetscCall(PetscSNPrintf(name, sizeof(name), "subdomain-%s-p%1.4" PetscInt_FMT ".vts", local_file_prefix, procid));
1616: for (II = 0; II < indent_level; II++) PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " "));
1618: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\" Source=\"%s\"/>\n", osx[i], oex[i] - 1, osy[j], oey[j] - 1, osz[k], oez[k] - 1, name));
1619: }
1620: }
1621: }
1622: PetscCall(PetscFree(olx));
1623: PetscCall(PetscFree(oly));
1624: PetscCall(PetscFree(olz));
1625: PetscCall(PetscFree(osx));
1626: PetscCall(PetscFree(osy));
1627: PetscCall(PetscFree(osz));
1628: PetscCall(PetscFree(oex));
1629: PetscCall(PetscFree(oey));
1630: PetscCall(PetscFree(oez));
1631: PetscFunctionReturn(PETSC_SUCCESS);
1632: }
1634: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da, const char file_prefix[], const char local_file_prefix[])
1635: {
1636: MPI_Comm comm;
1637: PetscMPIInt size, rank;
1638: char vtk_filename[PETSC_MAX_PATH_LEN];
1639: FILE *vtk_fp = NULL;
1640: const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
1641: PetscInt M, N, P, si, sj, sk, nx, ny, nz;
1642: PetscInt i, dofs;
1644: PetscFunctionBeginUser;
1645: /* only rank-0 generates this file */
1646: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
1647: PetscCallMPI(MPI_Comm_size(comm, &size));
1648: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1650: if (rank != 0) PetscFunctionReturn(PETSC_SUCCESS);
1652: /* create file name */
1653: PetscCall(PetscSNPrintf(vtk_filename, sizeof(vtk_filename), "%s.pvts", file_prefix));
1654: vtk_fp = fopen(vtk_filename, "w");
1655: PetscCheck(vtk_fp, PETSC_COMM_SELF, PETSC_ERR_SYS, "Cannot open file = %s ", vtk_filename);
1657: /* (VTK) generate pvts header */
1658: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<?xml version=\"1.0\"?>\n"));
1660: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n", byte_order));
1662: /* define size of the nodal mesh based on the cell DM */
1663: PetscCall(DMDAGetInfo(da, 0, &M, &N, &P, 0, 0, 0, &dofs, 0, 0, 0, 0, 0));
1664: PetscCall(DMDAGetGhostCorners(da, &si, &sj, &sk, &nx, &ny, &nz));
1665: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PStructuredGrid GhostLevel=\"1\" WholeExtent=\"%d %" PetscInt_FMT " %d %" PetscInt_FMT " %d %" PetscInt_FMT "\">\n", 0, M - 1, 0, N - 1, 0, P - 1)); /* note overlap = 1 for Q1 */
1667: /* DUMP THE CELL REFERENCES */
1668: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PCellData>\n"));
1669: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </PCellData>\n"));
1671: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PPoints>\n"));
1672: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PDataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"%d\"/>\n", NSD));
1673: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </PPoints>\n"));
1675: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PPointData>\n"));
1676: for (i = 0; i < dofs; i++) {
1677: const char *fieldname;
1678: PetscCall(DMDAGetFieldName(da, i, &fieldname));
1679: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n", fieldname));
1680: }
1681: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </PPointData>\n"));
1683: /* write out the parallel information */
1684: PetscCall(DAViewVTK_write_PieceExtend(vtk_fp, 2, da, local_file_prefix));
1686: /* close the file */
1687: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, " </PStructuredGrid>\n"));
1688: PetscCall(PetscFPrintf(PETSC_COMM_SELF, vtk_fp, "</VTKFile>\n"));
1690: if (vtk_fp) {
1691: fclose(vtk_fp);
1692: vtk_fp = NULL;
1693: }
1694: PetscFunctionReturn(PETSC_SUCCESS);
1695: }
1697: PetscErrorCode DAView3DPVTS(DM da, Vec x, const char NAME[])
1698: {
1699: char vts_filename[PETSC_MAX_PATH_LEN];
1700: char pvts_filename[PETSC_MAX_PATH_LEN];
1702: PetscFunctionBeginUser;
1703: PetscCall(PetscSNPrintf(vts_filename, sizeof(vts_filename), "%s-mesh", NAME));
1704: PetscCall(DAView_3DVTK_StructuredGrid_appended(da, x, vts_filename));
1706: PetscCall(PetscSNPrintf(pvts_filename, sizeof(pvts_filename), "%s-mesh", NAME));
1707: PetscCall(DAView_3DVTK_PStructuredGrid(da, pvts_filename, vts_filename));
1708: PetscFunctionReturn(PETSC_SUCCESS);
1709: }
1711: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp, PetscInt n, PetscReal rnorm, void *dummy)
1712: {
1713: PetscReal norms[4];
1714: Vec Br, v, w;
1715: Mat A;
1717: PetscFunctionBeginUser;
1718: PetscCall(KSPGetOperators(ksp, &A, NULL));
1719: PetscCall(MatCreateVecs(A, &w, &v));
1721: PetscCall(KSPBuildResidual(ksp, v, w, &Br));
1723: PetscCall(VecStrideNorm(Br, 0, NORM_2, &norms[0]));
1724: PetscCall(VecStrideNorm(Br, 1, NORM_2, &norms[1]));
1725: PetscCall(VecStrideNorm(Br, 2, NORM_2, &norms[2]));
1726: PetscCall(VecStrideNorm(Br, 3, NORM_2, &norms[3]));
1728: PetscCall(VecDestroy(&v));
1729: PetscCall(VecDestroy(&w));
1731: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%3" PetscInt_FMT " KSP Component U,V,W,P residual norm [ %1.12e, %1.12e, %1.12e, %1.12e ]\n", n, (double)norms[0], (double)norms[1], (double)norms[2], (double)norms[3]));
1732: PetscFunctionReturn(PETSC_SUCCESS);
1733: }
1735: static PetscErrorCode PCMGSetupViaCoarsen(PC pc, DM da_fine)
1736: {
1737: PetscInt nlevels, k;
1738: PETSC_UNUSED PetscInt finest;
1739: DM *da_list, *daclist;
1740: Mat R;
1742: PetscFunctionBeginUser;
1743: nlevels = 1;
1744: PetscCall(PetscOptionsGetInt(NULL, NULL, "-levels", &nlevels, 0));
1746: PetscCall(PetscMalloc1(nlevels, &da_list));
1747: for (k = 0; k < nlevels; k++) da_list[k] = NULL;
1748: PetscCall(PetscMalloc1(nlevels, &daclist));
1749: for (k = 0; k < nlevels; k++) daclist[k] = NULL;
1751: /* finest grid is nlevels - 1 */
1752: finest = nlevels - 1;
1753: daclist[0] = da_fine;
1754: PetscCall(PetscObjectReference((PetscObject)da_fine));
1755: PetscCall(DMCoarsenHierarchy(da_fine, nlevels - 1, &daclist[1]));
1756: for (k = 0; k < nlevels; k++) {
1757: da_list[k] = daclist[nlevels - 1 - k];
1758: PetscCall(DMDASetUniformCoordinates(da_list[k], 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
1759: }
1761: PetscCall(PCMGSetLevels(pc, nlevels, NULL));
1762: PetscCall(PCMGSetType(pc, PC_MG_MULTIPLICATIVE));
1763: PetscCall(PCMGSetGalerkin(pc, PC_MG_GALERKIN_PMAT));
1765: for (k = 1; k < nlevels; k++) {
1766: PetscCall(DMCreateInterpolation(da_list[k - 1], da_list[k], &R, NULL));
1767: PetscCall(PCMGSetInterpolation(pc, k, R));
1768: PetscCall(MatDestroy(&R));
1769: }
1771: /* tidy up */
1772: for (k = 0; k < nlevels; k++) PetscCall(DMDestroy(&da_list[k]));
1773: PetscCall(PetscFree(da_list));
1774: PetscCall(PetscFree(daclist));
1775: PetscFunctionReturn(PETSC_SUCCESS);
1776: }
1778: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx, PetscInt my, PetscInt mz)
1779: {
1780: DM da_Stokes;
1781: PetscInt u_dof, p_dof, dof, stencil_width;
1782: Mat A, B;
1783: DM vel_cda;
1784: Vec vel_coords;
1785: PetscInt p;
1786: Vec f, X, X1;
1787: DMDACoor3d ***_vel_coords;
1788: PetscInt its;
1789: KSP ksp_S;
1790: PetscInt model_definition = 0;
1791: PetscInt ei, ej, ek, sex, sey, sez, Imx, Jmy, Kmz;
1792: CellProperties cell_properties;
1793: PetscBool write_output = PETSC_FALSE, resolve = PETSC_FALSE;
1795: PetscFunctionBeginUser;
1796: PetscCall(PetscOptionsGetBool(NULL, NULL, "-resolve", &resolve, NULL));
1797: /* Generate the da for velocity and pressure */
1798: /* Num nodes in each direction is mx+1, my+1, mz+1 */
1799: u_dof = U_DOFS; /* Vx, Vy - velocities */
1800: p_dof = P_DOFS; /* p - pressure */
1801: dof = u_dof + p_dof;
1802: stencil_width = 1;
1803: PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, mz + 1, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, NULL, &da_Stokes));
1804: PetscCall(DMSetMatType(da_Stokes, MATAIJ));
1805: PetscCall(DMSetFromOptions(da_Stokes));
1806: PetscCall(DMSetUp(da_Stokes));
1807: PetscCall(DMDASetFieldName(da_Stokes, 0, "Vx"));
1808: PetscCall(DMDASetFieldName(da_Stokes, 1, "Vy"));
1809: PetscCall(DMDASetFieldName(da_Stokes, 2, "Vz"));
1810: PetscCall(DMDASetFieldName(da_Stokes, 3, "P"));
1812: /* unit box [0,1] x [0,1] x [0,1] */
1813: PetscCall(DMDASetUniformCoordinates(da_Stokes, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0));
1815: /* create quadrature point info for PDE coefficients */
1816: PetscCall(CellPropertiesCreate(da_Stokes, &cell_properties));
1818: /* interpolate the coordinates to quadrature points */
1819: PetscCall(DMGetCoordinateDM(da_Stokes, &vel_cda));
1820: PetscCall(DMGetCoordinatesLocal(da_Stokes, &vel_coords));
1821: PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));
1822: PetscCall(DMDAGetElementsCorners(da_Stokes, &sex, &sey, &sez));
1823: PetscCall(DMDAGetElementsSizes(da_Stokes, &Imx, &Jmy, &Kmz));
1824: for (ek = sez; ek < sez + Kmz; ek++) {
1825: for (ej = sey; ej < sey + Jmy; ej++) {
1826: for (ei = sex; ei < sex + Imx; ei++) {
1827: /* get coords for the element */
1828: PetscInt ngp;
1829: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
1830: PetscScalar el_coords[NSD * NODES_PER_EL];
1831: GaussPointCoefficients *cell;
1833: PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
1834: PetscCall(GetElementCoords3D(_vel_coords, ei, ej, ek, el_coords));
1835: ConstructGaussQuadrature3D(&ngp, gp_xi, gp_weight);
1837: for (p = 0; p < GAUSS_POINTS; p++) {
1838: PetscScalar xi_p[NSD], Ni_p[NODES_PER_EL];
1839: PetscScalar gp_x, gp_y, gp_z;
1840: PetscInt n;
1842: xi_p[0] = gp_xi[p][0];
1843: xi_p[1] = gp_xi[p][1];
1844: xi_p[2] = gp_xi[p][2];
1845: ShapeFunctionQ13D_Evaluate(xi_p, Ni_p);
1847: gp_x = gp_y = gp_z = 0.0;
1848: for (n = 0; n < NODES_PER_EL; n++) {
1849: gp_x = gp_x + Ni_p[n] * el_coords[NSD * n];
1850: gp_y = gp_y + Ni_p[n] * el_coords[NSD * n + 1];
1851: gp_z = gp_z + Ni_p[n] * el_coords[NSD * n + 2];
1852: }
1853: cell->gp_coords[NSD * p] = gp_x;
1854: cell->gp_coords[NSD * p + 1] = gp_y;
1855: cell->gp_coords[NSD * p + 2] = gp_z;
1856: }
1857: }
1858: }
1859: }
1861: PetscCall(PetscOptionsGetInt(NULL, NULL, "-model", &model_definition, NULL));
1863: switch (model_definition) {
1864: case 0: /* isoviscous */
1865: for (ek = sez; ek < sez + Kmz; ek++) {
1866: for (ej = sey; ej < sey + Jmy; ej++) {
1867: for (ei = sex; ei < sex + Imx; ei++) {
1868: GaussPointCoefficients *cell;
1870: PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
1871: for (p = 0; p < GAUSS_POINTS; p++) {
1872: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD * p]);
1873: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD * p + 1]);
1874: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD * p + 2]);
1876: cell->eta[p] = 1.0;
1878: cell->fx[p] = 0.0 * coord_x;
1879: cell->fy[p] = -PetscSinReal(2.2 * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1880: cell->fz[p] = 0.0 * coord_z;
1881: cell->hc[p] = 0.0;
1882: }
1883: }
1884: }
1885: }
1886: break;
1888: case 1: /* manufactured */
1889: for (ek = sez; ek < sez + Kmz; ek++) {
1890: for (ej = sey; ej < sey + Jmy; ej++) {
1891: for (ei = sex; ei < sex + Imx; ei++) {
1892: PetscReal eta, Fm[NSD], Fc, pos[NSD];
1893: GaussPointCoefficients *cell;
1895: PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
1896: for (p = 0; p < GAUSS_POINTS; p++) {
1897: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD * p]);
1898: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD * p + 1]);
1899: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD * p + 2]);
1901: pos[0] = coord_x;
1902: pos[1] = coord_y;
1903: pos[2] = coord_z;
1905: evaluate_MS_FrankKamentski(pos, NULL, NULL, &eta, Fm, &Fc);
1906: cell->eta[p] = eta;
1908: cell->fx[p] = Fm[0];
1909: cell->fy[p] = Fm[1];
1910: cell->fz[p] = Fm[2];
1911: cell->hc[p] = Fc;
1912: }
1913: }
1914: }
1915: }
1916: break;
1918: case 2: /* solcx */
1919: for (ek = sez; ek < sez + Kmz; ek++) {
1920: for (ej = sey; ej < sey + Jmy; ej++) {
1921: for (ei = sex; ei < sex + Imx; ei++) {
1922: GaussPointCoefficients *cell;
1924: PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
1925: for (p = 0; p < GAUSS_POINTS; p++) {
1926: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD * p]);
1927: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD * p + 1]);
1928: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD * p + 2]);
1930: cell->eta[p] = 1.0;
1932: cell->fx[p] = 0.0;
1933: cell->fy[p] = -PetscSinReal(3.0 * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1934: cell->fz[p] = 0.0 * coord_z;
1935: cell->hc[p] = 0.0;
1936: }
1937: }
1938: }
1939: }
1940: break;
1942: case 3: /* sinker */
1943: for (ek = sez; ek < sez + Kmz; ek++) {
1944: for (ej = sey; ej < sey + Jmy; ej++) {
1945: for (ei = sex; ei < sex + Imx; ei++) {
1946: GaussPointCoefficients *cell;
1948: PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
1949: for (p = 0; p < GAUSS_POINTS; p++) {
1950: PetscReal xp = PetscRealPart(cell->gp_coords[NSD * p]);
1951: PetscReal yp = PetscRealPart(cell->gp_coords[NSD * p + 1]);
1952: PetscReal zp = PetscRealPart(cell->gp_coords[NSD * p + 2]);
1954: cell->eta[p] = 1.0e-2;
1955: cell->fx[p] = 0.0;
1956: cell->fy[p] = 0.0;
1957: cell->fz[p] = 0.0;
1958: cell->hc[p] = 0.0;
1960: if ((PetscAbs(xp - 0.5) < 0.2) && (PetscAbs(yp - 0.5) < 0.2) && (PetscAbs(zp - 0.5) < 0.2)) {
1961: cell->eta[p] = 1.0;
1962: cell->fz[p] = 1.0;
1963: }
1964: }
1965: }
1966: }
1967: }
1968: break;
1970: case 4: /* subdomain jumps */
1971: for (ek = sez; ek < sez + Kmz; ek++) {
1972: for (ej = sey; ej < sey + Jmy; ej++) {
1973: for (ei = sex; ei < sex + Imx; ei++) {
1974: PetscReal opts_mag, opts_eta0;
1975: PetscInt px, py, pz;
1976: PetscBool jump;
1977: PetscMPIInt rr;
1978: GaussPointCoefficients *cell;
1980: opts_mag = 1.0;
1981: opts_eta0 = 1.e-2;
1983: PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_eta0", &opts_eta0, NULL));
1984: PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_magnitude", &opts_mag, NULL));
1985: PetscCall(DMDAGetInfo(da_Stokes, NULL, NULL, NULL, NULL, &px, &py, &pz, NULL, NULL, NULL, NULL, NULL, NULL));
1986: rr = PetscGlobalRank % (px * py);
1987: if (px % 2) jump = (PetscBool)(rr % 2);
1988: else jump = (PetscBool)((rr / px) % 2 ? rr % 2 : !(rr % 2));
1989: rr = PetscGlobalRank / (px * py);
1990: if (rr % 2) jump = (PetscBool)!jump;
1991: PetscCall(CellPropertiesGetCell(cell_properties, ei, ej, ek, &cell));
1992: for (p = 0; p < GAUSS_POINTS; p++) {
1993: PetscReal xp = PetscRealPart(cell->gp_coords[NSD * p]);
1994: PetscReal yp = PetscRealPart(cell->gp_coords[NSD * p + 1]);
1996: cell->eta[p] = jump ? PetscPowReal(10.0, opts_mag) : opts_eta0;
1997: cell->fx[p] = 0.0;
1998: cell->fy[p] = -PetscSinReal(2.2 * PETSC_PI * yp) * PetscCosReal(1.0 * PETSC_PI * xp);
1999: cell->fz[p] = 0.0;
2000: cell->hc[p] = 0.0;
2001: }
2002: }
2003: }
2004: }
2005: break;
2006: default:
2007: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "No default model is supported. Choose either -model {0,1,2,3}");
2008: }
2010: PetscCall(DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords));
2012: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
2013: PetscCall(DMCreateMatrix(da_Stokes, &A));
2014: PetscCall(DMCreateMatrix(da_Stokes, &B));
2015: PetscCall(DMCreateGlobalVector(da_Stokes, &X));
2016: PetscCall(DMCreateGlobalVector(da_Stokes, &f));
2018: /* assemble A11 */
2019: PetscCall(MatZeroEntries(A));
2020: PetscCall(MatZeroEntries(B));
2021: PetscCall(VecZeroEntries(f));
2023: PetscCall(AssembleA_Stokes(A, da_Stokes, cell_properties));
2024: PetscCall(MatViewFromOptions(A, NULL, "-amat_view"));
2025: PetscCall(AssembleA_PCStokes(B, da_Stokes, cell_properties));
2026: PetscCall(MatViewFromOptions(B, NULL, "-bmat_view"));
2027: /* build force vector */
2028: PetscCall(AssembleF_Stokes(f, da_Stokes, cell_properties));
2030: /* SOLVE */
2031: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_S));
2032: PetscCall(KSPSetOptionsPrefix(ksp_S, "stokes_"));
2033: PetscCall(KSPSetOperators(ksp_S, A, B));
2034: PetscCall(KSPSetFromOptions(ksp_S));
2036: {
2037: PC pc;
2038: const PetscInt ufields[] = {0, 1, 2}, pfields[] = {3};
2039: PetscCall(KSPGetPC(ksp_S, &pc));
2040: PetscCall(PCFieldSplitSetBlockSize(pc, 4));
2041: PetscCall(PCFieldSplitSetFields(pc, "u", 3, ufields, ufields));
2042: PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
2043: }
2045: {
2046: PC pc;
2047: PetscBool same = PETSC_FALSE;
2048: PetscCall(KSPGetPC(ksp_S, &pc));
2049: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCMG, &same));
2050: if (same) PetscCall(PCMGSetupViaCoarsen(pc, da_Stokes));
2051: }
2053: {
2054: PC pc;
2055: PetscBool same = PETSC_FALSE;
2056: PetscCall(KSPGetPC(ksp_S, &pc));
2057: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &same));
2058: if (same) PetscCall(KSPSetOperators(ksp_S, A, A));
2059: }
2061: {
2062: PetscBool stokes_monitor = PETSC_FALSE;
2063: PetscCall(PetscOptionsGetBool(NULL, NULL, "-stokes_ksp_monitor_blocks", &stokes_monitor, 0));
2064: if (stokes_monitor) PetscCall(KSPMonitorSet(ksp_S, KSPMonitorStokesBlocks, NULL, NULL));
2065: }
2067: if (resolve) {
2068: /* Test changing matrix data structure and resolve */
2069: PetscCall(VecDuplicate(X, &X1));
2070: }
2072: PetscCall(KSPSolve(ksp_S, f, X));
2073: if (resolve) {
2074: Mat C;
2075: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &C));
2076: PetscCall(KSPSetOperators(ksp_S, C, C));
2077: PetscCall(KSPSolve(ksp_S, f, X1));
2078: PetscCall(MatDestroy(&C));
2079: PetscCall(VecDestroy(&X1));
2080: }
2082: PetscCall(PetscOptionsGetBool(NULL, NULL, "-write_pvts", &write_output, NULL));
2083: if (write_output) PetscCall(DAView3DPVTS(da_Stokes, X, "up"));
2084: {
2085: PetscBool flg = PETSC_FALSE;
2086: char filename[PETSC_MAX_PATH_LEN];
2087: PetscCall(PetscOptionsGetString(NULL, NULL, "-write_binary", filename, sizeof(filename), &flg));
2088: if (flg) {
2089: PetscViewer viewer;
2090: /* PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer)); */
2091: PetscCall(PetscViewerVTKOpen(PETSC_COMM_WORLD, "ex42.vts", FILE_MODE_WRITE, &viewer));
2092: PetscCall(VecView(X, viewer));
2093: PetscCall(PetscViewerDestroy(&viewer));
2094: }
2095: }
2096: PetscCall(KSPGetIterationNumber(ksp_S, &its));
2098: /* verify */
2099: if (model_definition == 1) {
2100: DM da_Stokes_analytic;
2101: Vec X_analytic;
2103: PetscCall(DMDACreateManufacturedSolution(mx, my, mz, &da_Stokes_analytic, &X_analytic));
2104: if (write_output) PetscCall(DAView3DPVTS(da_Stokes_analytic, X_analytic, "ms"));
2105: PetscCall(DMDAIntegrateErrors3D(da_Stokes_analytic, X, X_analytic));
2106: if (write_output) PetscCall(DAView3DPVTS(da_Stokes, X, "up2"));
2107: PetscCall(DMDestroy(&da_Stokes_analytic));
2108: PetscCall(VecDestroy(&X_analytic));
2109: }
2111: PetscCall(KSPDestroy(&ksp_S));
2112: PetscCall(VecDestroy(&X));
2113: PetscCall(VecDestroy(&f));
2114: PetscCall(MatDestroy(&A));
2115: PetscCall(MatDestroy(&B));
2117: PetscCall(CellPropertiesDestroy(&cell_properties));
2118: PetscCall(DMDestroy(&da_Stokes));
2119: PetscFunctionReturn(PETSC_SUCCESS);
2120: }
2122: int main(int argc, char **args)
2123: {
2124: PetscInt mx, my, mz;
2126: PetscFunctionBeginUser;
2127: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
2128: mx = my = mz = 10;
2129: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
2130: my = mx;
2131: mz = mx;
2132: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
2133: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mz", &mz, NULL));
2134: PetscCall(solve_stokes_3d_coupled(mx, my, mz));
2135: PetscCall(PetscFinalize());
2136: return 0;
2137: }
2139: /*TEST
2141: test:
2142: args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type lu
2144: test:
2145: suffix: 2
2146: nsize: 3
2147: args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type redundant -stokes_redundant_pc_type lu
2149: test:
2150: suffix: bddc_stokes
2151: nsize: 6
2152: args: -mx 5 -my 4 -mz 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_dirichlet_pc_type svd -stokes_pc_bddc_neumann_pc_type svd -stokes_pc_bddc_coarse_redundant_pc_type svd
2154: test:
2155: suffix: bddc_stokes_deluxe
2156: nsize: 6
2157: args: -mx 5 -my 4 -mz 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_dirichlet_pc_type svd -stokes_pc_bddc_neumann_pc_type svd -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
2159: test:
2160: requires: !single
2161: suffix: bddc_stokes_subdomainjump_deluxe
2162: nsize: 9
2163: args: -model 4 -jump_magnitude 4 -mx 6 -my 6 -mz 2 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -stokes_pc_bddc_use_deluxe_scaling -stokes_sub_schurs_posdef 0 -stokes_sub_schurs_symmetric -stokes_sub_schurs_mat_solver_type petsc -stokes_pc_bddc_schur_layers 1
2165: test:
2166: requires: !single
2167: suffix: 3
2168: args: -stokes_ksp_converged_reason -stokes_pc_type fieldsplit -resolve
2170: test:
2171: suffix: tut_1
2172: nsize: 4
2173: requires: !single
2174: args: -stokes_ksp_monitor
2176: test:
2177: suffix: tut_2
2178: nsize: 4
2179: requires: !single
2180: args: -stokes_ksp_monitor -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type schur
2182: test:
2183: suffix: tut_3
2184: nsize: 4
2185: requires: !single
2186: args: -mx 20 -stokes_ksp_monitor -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type schur
2188: TEST*/