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;
56: PetscNew(&cells);
58: DMDAGetElementsCorners(da_stokes,&sex,&sey,&sez);
59: 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: PetscMalloc1(mx*my*mz,&cells->gpc);
70: *C = cells;
71: return 0;
72: }
74: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
75: {
76: CellProperties cells;
79: if (!C) return 0;
80: cells = *C;
81: PetscFree(cells->gpc);
82: PetscFree(cells);
83: *C = NULL;
84: return 0;
85: }
87: PetscErrorCode CellPropertiesGetCell(CellProperties C,PetscInt II,PetscInt J,PetscInt K,GaussPointCoefficients **G)
88: {
90: *G = &C->gpc[(II-C->sex) + (J-C->sey)*C->mx + (K-C->sez)*C->mx*C->my];
91: return 0;
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; JJ[0][1] = J01; JJ[0][2] = J02;
208: JJ[1][0] = J10; JJ[1][1] = J11; JJ[1][2] = J12;
209: JJ[2][0] = J20; JJ[2][1] = J21; JJ[2][2] = J22;
211: matrix_inverse_3x3(JJ,iJ);
213: *det_J = J00*J11*J22 - J00*J12*J21 - J10*J01*J22 + J10*J02*J21 + J20*J01*J12 - J20*J02*J11;
215: for (n=0; n<NODES_PER_EL; n++) {
216: GNx[0][n] = GNi[0][n]*iJ[0][0] + GNi[1][n]*iJ[0][1] + GNi[2][n]*iJ[0][2];
217: GNx[1][n] = GNi[0][n]*iJ[1][0] + GNi[1][n]*iJ[1][1] + GNi[2][n]*iJ[1][2];
218: GNx[2][n] = GNi[0][n]*iJ[2][0] + GNi[1][n]*iJ[2][1] + GNi[2][n]*iJ[2][2];
219: }
220: }
222: static void ConstructGaussQuadrature3D(PetscInt *ngp,PetscScalar gp_xi[][NSD],PetscScalar gp_weight[])
223: {
224: *ngp = 8;
225: gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919; gp_xi[0][2] = -0.57735026919;
226: gp_xi[1][0] = -0.57735026919; gp_xi[1][1] = 0.57735026919; gp_xi[1][2] = -0.57735026919;
227: gp_xi[2][0] = 0.57735026919; gp_xi[2][1] = 0.57735026919; gp_xi[2][2] = -0.57735026919;
228: gp_xi[3][0] = 0.57735026919; gp_xi[3][1] = -0.57735026919; gp_xi[3][2] = -0.57735026919;
230: gp_xi[4][0] = -0.57735026919; gp_xi[4][1] = -0.57735026919; gp_xi[4][2] = 0.57735026919;
231: gp_xi[5][0] = -0.57735026919; gp_xi[5][1] = 0.57735026919; gp_xi[5][2] = 0.57735026919;
232: gp_xi[6][0] = 0.57735026919; gp_xi[6][1] = 0.57735026919; gp_xi[6][2] = 0.57735026919;
233: gp_xi[7][0] = 0.57735026919; gp_xi[7][1] = -0.57735026919; gp_xi[7][2] = 0.57735026919;
235: gp_weight[0] = 1.0;
236: gp_weight[1] = 1.0;
237: gp_weight[2] = 1.0;
238: gp_weight[3] = 1.0;
240: gp_weight[4] = 1.0;
241: gp_weight[5] = 1.0;
242: gp_weight[6] = 1.0;
243: gp_weight[7] = 1.0;
244: }
246: /*
247: i,j are the element indices
248: The unknown is a vector quantity.
249: The s[].c is used to indicate the degree of freedom.
250: */
251: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j,PetscInt k)
252: {
253: PetscInt n;
256: /* velocity */
257: n = 0;
258: /* node 0 */
259: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++; /* Vx0 */
260: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++; /* Vy0 */
261: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++; /* Vz0 */
263: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
264: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
265: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
267: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
268: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
269: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
271: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++;
272: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++;
273: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++;
275: /* */
276: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++; /* Vx4 */
277: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++; /* Vy4 */
278: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++; /* Vz4 */
280: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
281: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
282: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
284: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
285: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
286: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
288: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++;
289: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++;
290: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++;
292: /* pressure */
293: n = 0;
295: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++; /* P0 */
296: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
297: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
298: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++;
300: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++; /* P0 */
301: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
302: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
303: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++;
304: return 0;
305: }
307: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords,PetscInt i,PetscInt j,PetscInt k,PetscScalar el_coord[])
308: {
310: /* get coords for the element */
311: el_coord[0] = coords[k][j][i].x;
312: el_coord[1] = coords[k][j][i].y;
313: el_coord[2] = coords[k][j][i].z;
315: el_coord[3] = coords[k][j+1][i].x;
316: el_coord[4] = coords[k][j+1][i].y;
317: el_coord[5] = coords[k][j+1][i].z;
319: el_coord[6] = coords[k][j+1][i+1].x;
320: el_coord[7] = coords[k][j+1][i+1].y;
321: el_coord[8] = coords[k][j+1][i+1].z;
323: el_coord[9] = coords[k][j][i+1].x;
324: el_coord[10] = coords[k][j][i+1].y;
325: el_coord[11] = coords[k][j][i+1].z;
327: el_coord[12] = coords[k+1][j][i].x;
328: el_coord[13] = coords[k+1][j][i].y;
329: el_coord[14] = coords[k+1][j][i].z;
331: el_coord[15] = coords[k+1][j+1][i].x;
332: el_coord[16] = coords[k+1][j+1][i].y;
333: el_coord[17] = coords[k+1][j+1][i].z;
335: el_coord[18] = coords[k+1][j+1][i+1].x;
336: el_coord[19] = coords[k+1][j+1][i+1].y;
337: el_coord[20] = coords[k+1][j+1][i+1].z;
339: el_coord[21] = coords[k+1][j][i+1].x;
340: el_coord[22] = coords[k+1][j][i+1].y;
341: el_coord[23] = coords[k+1][j][i+1].z;
342: return 0;
343: }
345: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field,PetscInt i,PetscInt j,PetscInt k,StokesDOF nodal_fields[])
346: {
348: /* get the nodal fields for u */
349: nodal_fields[0].u_dof = field[k][j][i].u_dof;
350: nodal_fields[0].v_dof = field[k][j][i].v_dof;
351: nodal_fields[0].w_dof = field[k][j][i].w_dof;
353: nodal_fields[1].u_dof = field[k][j+1][i].u_dof;
354: nodal_fields[1].v_dof = field[k][j+1][i].v_dof;
355: nodal_fields[1].w_dof = field[k][j+1][i].w_dof;
357: nodal_fields[2].u_dof = field[k][j+1][i+1].u_dof;
358: nodal_fields[2].v_dof = field[k][j+1][i+1].v_dof;
359: nodal_fields[2].w_dof = field[k][j+1][i+1].w_dof;
361: nodal_fields[3].u_dof = field[k][j][i+1].u_dof;
362: nodal_fields[3].v_dof = field[k][j][i+1].v_dof;
363: nodal_fields[3].w_dof = field[k][j][i+1].w_dof;
365: nodal_fields[4].u_dof = field[k+1][j][i].u_dof;
366: nodal_fields[4].v_dof = field[k+1][j][i].v_dof;
367: nodal_fields[4].w_dof = field[k+1][j][i].w_dof;
369: nodal_fields[5].u_dof = field[k+1][j+1][i].u_dof;
370: nodal_fields[5].v_dof = field[k+1][j+1][i].v_dof;
371: nodal_fields[5].w_dof = field[k+1][j+1][i].w_dof;
373: nodal_fields[6].u_dof = field[k+1][j+1][i+1].u_dof;
374: nodal_fields[6].v_dof = field[k+1][j+1][i+1].v_dof;
375: nodal_fields[6].w_dof = field[k+1][j+1][i+1].w_dof;
377: nodal_fields[7].u_dof = field[k+1][j][i+1].u_dof;
378: nodal_fields[7].v_dof = field[k+1][j][i+1].v_dof;
379: nodal_fields[7].w_dof = field[k+1][j][i+1].w_dof;
381: /* get the nodal fields for p */
382: nodal_fields[0].p_dof = field[k][j][i].p_dof;
383: nodal_fields[1].p_dof = field[k][j+1][i].p_dof;
384: nodal_fields[2].p_dof = field[k][j+1][i+1].p_dof;
385: nodal_fields[3].p_dof = field[k][j][i+1].p_dof;
387: nodal_fields[4].p_dof = field[k+1][j][i].p_dof;
388: nodal_fields[5].p_dof = field[k+1][j+1][i].p_dof;
389: nodal_fields[6].p_dof = field[k+1][j+1][i+1].p_dof;
390: nodal_fields[7].p_dof = field[k+1][j][i+1].p_dof;
391: return 0;
392: }
394: 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)
395: {
396: PetscInt ij;
397: PETSC_UNUSED PetscInt r,c,nr,nc;
399: nr = w_NPE*w_dof;
400: nc = u_NPE*u_dof;
402: r = w_dof*wi+wd;
403: c = u_dof*ui+ud;
405: ij = r*nc+c;
407: return ij;
408: }
410: static PetscErrorCode DMDASetValuesLocalStencil3D_ADD_VALUES(StokesDOF ***fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
411: {
412: PetscInt n,II,J,K;
415: for (n = 0; n<NODES_PER_EL; n++) {
416: II = u_eqn[NSD*n].i;
417: J = u_eqn[NSD*n].j;
418: K = u_eqn[NSD*n].k;
420: fields_F[K][J][II].u_dof = fields_F[K][J][II].u_dof+Fe_u[NSD*n];
422: II = u_eqn[NSD*n+1].i;
423: J = u_eqn[NSD*n+1].j;
424: K = u_eqn[NSD*n+1].k;
426: fields_F[K][J][II].v_dof = fields_F[K][J][II].v_dof+Fe_u[NSD*n+1];
428: II = u_eqn[NSD*n+2].i;
429: J = u_eqn[NSD*n+2].j;
430: K = u_eqn[NSD*n+2].k;
431: fields_F[K][J][II].w_dof = fields_F[K][J][II].w_dof+Fe_u[NSD*n+2];
433: II = p_eqn[n].i;
434: J = p_eqn[n].j;
435: K = p_eqn[n].k;
437: fields_F[K][J][II].p_dof = fields_F[K][J][II].p_dof+Fe_p[n];
439: }
440: return 0;
441: }
443: static void FormStressOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
444: {
445: PetscInt ngp;
446: PetscScalar gp_xi[GAUSS_POINTS][NSD];
447: PetscScalar gp_weight[GAUSS_POINTS];
448: PetscInt p,i,j,k;
449: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
450: PetscScalar J_p,tildeD[6];
451: PetscScalar B[6][U_DOFS*NODES_PER_EL];
452: const PetscInt nvdof = U_DOFS*NODES_PER_EL;
454: /* define quadrature rule */
455: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
457: /* evaluate integral */
458: for (p = 0; p < ngp; p++) {
459: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
460: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
462: for (i = 0; i < NODES_PER_EL; i++) {
463: PetscScalar d_dx_i = GNx_p[0][i];
464: PetscScalar d_dy_i = GNx_p[1][i];
465: PetscScalar d_dz_i = GNx_p[2][i];
467: B[0][3*i] = d_dx_i; B[0][3*i+1] = 0.0; B[0][3*i+2] = 0.0;
468: B[1][3*i] = 0.0; B[1][3*i+1] = d_dy_i; B[1][3*i+2] = 0.0;
469: B[2][3*i] = 0.0; B[2][3*i+1] = 0.0; B[2][3*i+2] = d_dz_i;
471: B[3][3*i] = d_dy_i; B[3][3*i+1] = d_dx_i; B[3][3*i+2] = 0.0; /* e_xy */
472: B[4][3*i] = d_dz_i; B[4][3*i+1] = 0.0; B[4][3*i+2] = d_dx_i; /* e_xz */
473: B[5][3*i] = 0.0; B[5][3*i+1] = d_dz_i; B[5][3*i+2] = d_dy_i; /* e_yz */
474: }
476: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
477: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
478: tildeD[2] = 2.0*gp_weight[p]*J_p*eta[p];
480: tildeD[3] = gp_weight[p]*J_p*eta[p];
481: tildeD[4] = gp_weight[p]*J_p*eta[p];
482: tildeD[5] = gp_weight[p]*J_p*eta[p];
484: /* form Bt tildeD B */
485: /*
486: Ke_ij = Bt_ik . D_kl . B_lj
487: = B_ki . D_kl . B_lj
488: = B_ki . D_kk . B_kj
489: */
490: for (i = 0; i < nvdof; i++) {
491: for (j = i; j < nvdof; j++) {
492: for (k = 0; k < 6; k++) {
493: Ke[i*nvdof+j] += B[k][i]*tildeD[k]*B[k][j];
494: }
495: }
496: }
498: }
499: /* fill lower triangular part */
500: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
501: for (i = 0; i < nvdof; i++) {
502: for (j = i; j < nvdof; j++) {
503: Ke[j*nvdof+i] = Ke[i*nvdof+j];
504: }
505: }
506: #endif
507: }
509: static void FormGradientOperatorQ13D(PetscScalar Ke[],PetscScalar coords[])
510: {
511: PetscInt ngp;
512: PetscScalar gp_xi[GAUSS_POINTS][NSD];
513: PetscScalar gp_weight[GAUSS_POINTS];
514: PetscInt p,i,j,di;
515: PetscScalar Ni_p[NODES_PER_EL];
516: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
517: PetscScalar J_p,fac;
519: /* define quadrature rule */
520: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
522: /* evaluate integral */
523: for (p = 0; p < ngp; p++) {
524: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
525: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
526: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
527: fac = gp_weight[p]*J_p;
529: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
530: for (di = 0; di < NSD; di++) { /* u dofs */
531: for (j = 0; j < NODES_PER_EL; j++) { /* p nodes, p dofs = 1 (ie no loop) */
532: PetscInt IJ;
533: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,3,j,0,NODES_PER_EL,1);
535: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
536: }
537: }
538: }
539: }
540: }
542: static void FormDivergenceOperatorQ13D(PetscScalar De[],PetscScalar coords[])
543: {
544: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
545: PetscInt i,j;
546: PetscInt nr_g,nc_g;
548: PetscMemzero(Ge,sizeof(Ge));
549: FormGradientOperatorQ13D(Ge,coords);
551: nr_g = U_DOFS*NODES_PER_EL;
552: nc_g = P_DOFS*NODES_PER_EL;
554: for (i = 0; i < nr_g; i++) {
555: for (j = 0; j < nc_g; j++) {
556: De[nr_g*j+i] = Ge[nc_g*i+j];
557: }
558: }
559: }
561: static void FormStabilisationOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
562: {
563: PetscInt ngp;
564: PetscScalar gp_xi[GAUSS_POINTS][NSD];
565: PetscScalar gp_weight[GAUSS_POINTS];
566: PetscInt p,i,j;
567: PetscScalar Ni_p[NODES_PER_EL];
568: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
569: PetscScalar J_p,fac,eta_avg;
571: /* define quadrature rule */
572: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
574: /* evaluate integral */
575: for (p = 0; p < ngp; p++) {
576: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
577: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
578: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
579: fac = gp_weight[p]*J_p;
580: /*
581: for (i = 0; i < NODES_PER_EL; i++) {
582: for (j = i; j < NODES_PER_EL; j++) {
583: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
584: }
585: }
586: */
588: for (i = 0; i < NODES_PER_EL; i++) {
589: for (j = 0; j < NODES_PER_EL; j++) {
590: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
591: }
592: }
593: }
595: /* scale */
596: eta_avg = 0.0;
597: for (p = 0; p < ngp; p++) eta_avg += eta[p];
598: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
599: fac = 1.0/eta_avg;
601: /*
602: for (i = 0; i < NODES_PER_EL; i++) {
603: for (j = i; j < NODES_PER_EL; j++) {
604: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
605: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
606: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
607: #endif
608: }
609: }
610: */
612: for (i = 0; i < NODES_PER_EL; i++) {
613: for (j = 0; j < NODES_PER_EL; j++) {
614: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
615: }
616: }
617: }
619: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
620: {
621: PetscInt ngp;
622: PetscScalar gp_xi[GAUSS_POINTS][NSD];
623: PetscScalar gp_weight[GAUSS_POINTS];
624: PetscInt p,i,j;
625: PetscScalar Ni_p[NODES_PER_EL];
626: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
627: PetscScalar J_p,fac,eta_avg;
629: /* define quadrature rule */
630: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
632: /* evaluate integral */
633: for (p = 0; p < ngp; p++) {
634: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
635: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
636: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
637: fac = gp_weight[p]*J_p;
639: /*
640: for (i = 0; i < NODES_PER_EL; i++) {
641: for (j = i; j < NODES_PER_EL; j++) {
642: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
643: }
644: }
645: */
647: for (i = 0; i < NODES_PER_EL; i++) {
648: for (j = 0; j < NODES_PER_EL; j++) {
649: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
650: }
651: }
652: }
654: /* scale */
655: eta_avg = 0.0;
656: for (p = 0; p < ngp; p++) eta_avg += eta[p];
657: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
658: fac = 1.0/eta_avg;
659: /*
660: for (i = 0; i < NODES_PER_EL; i++) {
661: for (j = i; j < NODES_PER_EL; j++) {
662: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
663: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
664: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
665: #endif
666: }
667: }
668: */
670: for (i = 0; i < NODES_PER_EL; i++) {
671: for (j = 0; j < NODES_PER_EL; j++) {
672: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
673: }
674: }
675: }
677: static void FormMomentumRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[],PetscScalar fz[])
678: {
679: PetscInt ngp;
680: PetscScalar gp_xi[GAUSS_POINTS][NSD];
681: PetscScalar gp_weight[GAUSS_POINTS];
682: PetscInt p,i;
683: PetscScalar Ni_p[NODES_PER_EL];
684: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
685: PetscScalar J_p,fac;
687: /* define quadrature rule */
688: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
690: /* evaluate integral */
691: for (p = 0; p < ngp; p++) {
692: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
693: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
694: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
695: fac = gp_weight[p]*J_p;
697: for (i = 0; i < NODES_PER_EL; i++) {
698: Fe[NSD*i] -= fac*Ni_p[i]*fx[p];
699: Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
700: Fe[NSD*i+2] -= fac*Ni_p[i]*fz[p];
701: }
702: }
703: }
705: static void FormContinuityRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar hc[])
706: {
707: PetscInt ngp;
708: PetscScalar gp_xi[GAUSS_POINTS][NSD];
709: PetscScalar gp_weight[GAUSS_POINTS];
710: PetscInt p,i;
711: PetscScalar Ni_p[NODES_PER_EL];
712: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
713: PetscScalar J_p,fac;
715: /* define quadrature rule */
716: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
718: /* evaluate integral */
719: for (p = 0; p < ngp; p++) {
720: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
721: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
722: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
723: fac = gp_weight[p]*J_p;
725: for (i = 0; i < NODES_PER_EL; i++) Fe[i] -= fac*Ni_p[i]*hc[p];
726: }
727: }
729: #define _ZERO_ROWCOL_i(A,i) { \
730: PetscInt KK; \
731: PetscScalar tmp = A[24*(i)+(i)]; \
732: for (KK=0;KK<24;KK++) A[24*(i)+KK]=0.0; \
733: for (KK=0;KK<24;KK++) A[24*KK+(i)]=0.0; \
734: A[24*(i)+(i)] = tmp;} \
736: #define _ZERO_ROW_i(A,i) { \
737: PetscInt KK; \
738: for (KK=0;KK<8;KK++) A[8*(i)+KK]=0.0;}
740: #define _ZERO_COL_i(A,i) { \
741: PetscInt KK; \
742: for (KK=0;KK<8;KK++) A[24*KK+(i)]=0.0;}
744: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,CellProperties cell_properties)
745: {
746: DM cda;
747: Vec coords;
748: DMDACoor3d ***_coords;
749: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
750: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
751: PetscInt sex,sey,sez,mx,my,mz;
752: PetscInt ei,ej,ek;
753: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
754: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
755: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
756: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
757: PetscScalar el_coords[NODES_PER_EL*NSD];
758: GaussPointCoefficients *props;
759: PetscScalar *prop_eta;
760: PetscInt n,M,N,P;
763: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
764: /* setup for coords */
765: DMGetCoordinateDM(stokes_da,&cda);
766: DMGetCoordinatesLocal(stokes_da,&coords);
767: DMDAVecGetArray(cda,coords,&_coords);
768: DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
769: DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
770: for (ek = sez; ek < sez+mz; ek++) {
771: for (ej = sey; ej < sey+my; ej++) {
772: for (ei = sex; ei < sex+mx; ei++) {
773: /* get coords for the element */
774: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
775: /* get cell properties */
776: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
777: /* get coefficients for the element */
778: prop_eta = props->eta;
780: /* initialise element stiffness matrix */
781: PetscMemzero(Ae,sizeof(Ae));
782: PetscMemzero(Ge,sizeof(Ge));
783: PetscMemzero(De,sizeof(De));
784: PetscMemzero(Ce,sizeof(Ce));
786: /* form element stiffness matrix */
787: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
788: FormGradientOperatorQ13D(Ge,el_coords);
789: /*#if defined(ASSEMBLE_LOWER_TRIANGULAR)*/
790: FormDivergenceOperatorQ13D(De,el_coords);
791: /*#endif*/
792: FormStabilisationOperatorQ13D(Ce,el_coords,prop_eta);
794: /* insert element matrix into global matrix */
795: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
797: for (n=0; n<NODES_PER_EL; n++) {
798: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
799: _ZERO_ROWCOL_i(Ae,3*n);
800: _ZERO_ROW_i (Ge,3*n);
801: _ZERO_COL_i (De,3*n);
802: }
804: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
805: _ZERO_ROWCOL_i(Ae,3*n+1);
806: _ZERO_ROW_i (Ge,3*n+1);
807: _ZERO_COL_i (De,3*n+1);
808: }
810: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
811: _ZERO_ROWCOL_i(Ae,3*n+2);
812: _ZERO_ROW_i (Ge,3*n+2);
813: _ZERO_COL_i (De,3*n+2);
814: }
815: }
816: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
817: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
818: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
819: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
820: }
821: }
822: }
823: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
824: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
826: DMDAVecRestoreArray(cda,coords,&_coords);
828: return 0;
829: }
831: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,CellProperties cell_properties)
832: {
833: DM cda;
834: Vec coords;
835: DMDACoor3d ***_coords;
836: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
837: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
838: PetscInt sex,sey,sez,mx,my,mz;
839: PetscInt ei,ej,ek;
840: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
841: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
842: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
843: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
844: PetscScalar el_coords[NODES_PER_EL*NSD];
845: GaussPointCoefficients *props;
846: PetscScalar *prop_eta;
847: PetscInt n,M,N,P;
850: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
851: /* setup for coords */
852: DMGetCoordinateDM(stokes_da,&cda);
853: DMGetCoordinatesLocal(stokes_da,&coords);
854: DMDAVecGetArray(cda,coords,&_coords);
856: DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
857: DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
858: for (ek = sez; ek < sez+mz; ek++) {
859: for (ej = sey; ej < sey+my; ej++) {
860: for (ei = sex; ei < sex+mx; ei++) {
861: /* get coords for the element */
862: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
863: /* get cell properties */
864: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
865: /* get coefficients for the element */
866: prop_eta = props->eta;
868: /* initialise element stiffness matrix */
869: PetscMemzero(Ae,sizeof(Ae));
870: PetscMemzero(Ge,sizeof(Ge));
871: PetscMemzero(De,sizeof(De));
872: PetscMemzero(Ce,sizeof(Ce));
874: /* form element stiffness matrix */
875: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
876: FormGradientOperatorQ13D(Ge,el_coords);
877: /* FormDivergenceOperatorQ13D(De,el_coords); */
878: FormScaledMassMatrixOperatorQ13D(Ce,el_coords,prop_eta);
880: /* insert element matrix into global matrix */
881: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
883: for (n=0; n<NODES_PER_EL; n++) {
884: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
885: _ZERO_ROWCOL_i(Ae,3*n);
886: _ZERO_ROW_i (Ge,3*n);
887: _ZERO_COL_i (De,3*n);
888: }
890: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
891: _ZERO_ROWCOL_i(Ae,3*n+1);
892: _ZERO_ROW_i (Ge,3*n+1);
893: _ZERO_COL_i (De,3*n+1);
894: }
896: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
897: _ZERO_ROWCOL_i(Ae,3*n+2);
898: _ZERO_ROW_i (Ge,3*n+2);
899: _ZERO_COL_i (De,3*n+2);
900: }
901: }
902: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
903: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
904: /*MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);*/
905: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
906: }
907: }
908: }
909: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
910: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
912: DMDAVecRestoreArray(cda,coords,&_coords);
913: return 0;
914: }
916: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,CellProperties cell_properties)
917: {
918: DM cda;
919: Vec coords;
920: DMDACoor3d ***_coords;
921: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
922: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
923: PetscInt sex,sey,sez,mx,my,mz;
924: PetscInt ei,ej,ek;
925: PetscScalar Fe[NODES_PER_EL*U_DOFS];
926: PetscScalar He[NODES_PER_EL*P_DOFS];
927: PetscScalar el_coords[NODES_PER_EL*NSD];
928: GaussPointCoefficients *props;
929: PetscScalar *prop_fx,*prop_fy,*prop_fz,*prop_hc;
930: Vec local_F;
931: StokesDOF ***ff;
932: PetscInt n,M,N,P;
935: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
936: /* setup for coords */
937: DMGetCoordinateDM(stokes_da,&cda);
938: DMGetCoordinatesLocal(stokes_da,&coords);
939: DMDAVecGetArray(cda,coords,&_coords);
941: /* get access to the vector */
942: DMGetLocalVector(stokes_da,&local_F);
943: VecZeroEntries(local_F);
944: DMDAVecGetArray(stokes_da,local_F,&ff);
945: DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
946: DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
947: for (ek = sez; ek < sez+mz; ek++) {
948: for (ej = sey; ej < sey+my; ej++) {
949: for (ei = sex; ei < sex+mx; ei++) {
950: /* get coords for the element */
951: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
952: /* get cell properties */
953: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
954: /* get coefficients for the element */
955: prop_fx = props->fx;
956: prop_fy = props->fy;
957: prop_fz = props->fz;
958: prop_hc = props->hc;
960: /* initialise element stiffness matrix */
961: PetscMemzero(Fe,sizeof(Fe));
962: PetscMemzero(He,sizeof(He));
964: /* form element stiffness matrix */
965: FormMomentumRhsQ13D(Fe,el_coords,prop_fx,prop_fy,prop_fz);
966: FormContinuityRhsQ13D(He,el_coords,prop_hc);
968: /* insert element matrix into global matrix */
969: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
971: for (n=0; n<NODES_PER_EL; n++) {
972: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) Fe[3*n] = 0.0;
974: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) Fe[3*n+1] = 0.0;
976: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) Fe[3*n+2] = 0.0;
977: }
979: DMDASetValuesLocalStencil3D_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
980: }
981: }
982: }
983: DMDAVecRestoreArray(stokes_da,local_F,&ff);
984: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
985: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
986: DMRestoreLocalVector(stokes_da,&local_F);
988: DMDAVecRestoreArray(cda,coords,&_coords);
989: return 0;
990: }
992: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta,PetscReal *MX,PetscReal *MY,PetscReal *MZ)
993: {
994: *theta = 0.0;
995: *MX = 2.0 * PETSC_PI;
996: *MY = 2.0 * PETSC_PI;
997: *MZ = 2.0 * PETSC_PI;
998: }
999: static void evaluate_MS_FrankKamentski(PetscReal pos[],PetscReal v[],PetscReal *p,PetscReal *eta,PetscReal Fm[],PetscReal *Fc)
1000: {
1001: PetscReal x,y,z;
1002: PetscReal theta,MX,MY,MZ;
1004: evaluate_MS_FrankKamentski_constants(&theta,&MX,&MY,&MZ);
1005: x = pos[0];
1006: y = pos[1];
1007: z = pos[2];
1008: if (v) {
1009: /*
1010: v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1011: v[1] = z*cos(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1012: v[2] = -(PetscPowRealInt(x,3) + PetscPowRealInt(y,3))*PetscSinReal(2.0*PETSC_PI*z);
1013: */
1014: /*
1015: v[0] = PetscSinReal(PETSC_PI*x);
1016: v[1] = PetscSinReal(PETSC_PI*y);
1017: v[2] = PetscSinReal(PETSC_PI*z);
1018: */
1019: v[0] = PetscPowRealInt(z,3)*PetscExpReal(y)*PetscSinReal(PETSC_PI*x);
1020: v[1] = z*PetscCosReal(2.0*PETSC_PI*x)*PetscExpReal(-y)*PetscSinReal(PETSC_PI*y);
1021: 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;
1022: }
1023: if (p) *p = PetscPowRealInt(x,2) + PetscPowRealInt(y,2) + PetscPowRealInt(z,2);
1024: if (eta) {
1025: /* eta = PetscExpReal(-theta*(1.0 - y - 0.1*PetscCosReal(MX*x)*PetscCosReal(MZ*z)*PetscSinReal(MY*y)));*/
1026: *eta = 1.0;
1027: }
1028: if (Fm) {
1029: /*
1030: 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))) ;
1031: 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)));
1032: 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))) ;
1033: */
1034: /*
1035: Fm[0]=-2*x - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*x);
1036: Fm[1]=-2*y - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*y);
1037: Fm[2]=-2*z - 2.0*PetscPowRealInt(PETSC_PI,2)*PetscSinReal(PETSC_PI*z);
1038: */
1039: /*
1040: 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) ;
1041: 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);
1042: 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) ;
1043: */
1044: 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);
1045: 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);
1046: 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);
1047: }
1048: if (Fc) {
1049: /* 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) ;*/
1050: /* Fc = PETSC_PI*PetscCosReal(PETSC_PI*x) + PETSC_PI*PetscCosReal(PETSC_PI*y) + PETSC_PI*PetscCosReal(PETSC_PI*z);*/
1051: /* 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);*/
1052: *Fc = 0.0;
1053: }
1054: }
1056: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx,PetscInt my,PetscInt mz,DM *_da,Vec *_X)
1057: {
1058: DM da,cda;
1059: Vec X;
1060: StokesDOF ***_stokes;
1061: Vec coords;
1062: DMDACoor3d ***_coords;
1063: PetscInt si,sj,sk,ei,ej,ek,i,j,k;
1066: PetscCall(DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1067: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,4,1,NULL,NULL,NULL,&da));
1068: DMSetFromOptions(da);
1069: DMSetUp(da);
1070: DMDASetFieldName(da,0,"anlytic_Vx");
1071: DMDASetFieldName(da,1,"anlytic_Vy");
1072: DMDASetFieldName(da,2,"anlytic_Vz");
1073: DMDASetFieldName(da,3,"analytic_P");
1075: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
1077: DMGetCoordinatesLocal(da,&coords);
1078: DMGetCoordinateDM(da,&cda);
1079: DMDAVecGetArray(cda,coords,&_coords);
1081: DMCreateGlobalVector(da,&X);
1082: DMDAVecGetArray(da,X,&_stokes);
1084: DMDAGetCorners(da,&si,&sj,&sk,&ei,&ej,&ek);
1085: for (k = sk; k < sk+ek; k++) {
1086: for (j = sj; j < sj+ej; j++) {
1087: for (i = si; i < si+ei; i++) {
1088: PetscReal pos[NSD],pressure,vel[NSD];
1090: pos[0] = PetscRealPart(_coords[k][j][i].x);
1091: pos[1] = PetscRealPart(_coords[k][j][i].y);
1092: pos[2] = PetscRealPart(_coords[k][j][i].z);
1094: evaluate_MS_FrankKamentski(pos,vel,&pressure,NULL,NULL,NULL);
1096: _stokes[k][j][i].u_dof = vel[0];
1097: _stokes[k][j][i].v_dof = vel[1];
1098: _stokes[k][j][i].w_dof = vel[2];
1099: _stokes[k][j][i].p_dof = pressure;
1100: }
1101: }
1102: }
1103: DMDAVecRestoreArray(da,X,&_stokes);
1104: DMDAVecRestoreArray(cda,coords,&_coords);
1106: *_da = da;
1107: *_X = X;
1108: return 0;
1109: }
1111: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da,Vec X,Vec X_analytic)
1112: {
1113: DM cda;
1114: Vec coords,X_analytic_local,X_local;
1115: DMDACoor3d ***_coords;
1116: PetscInt sex,sey,sez,mx,my,mz;
1117: PetscInt ei,ej,ek;
1118: PetscScalar el_coords[NODES_PER_EL*NSD];
1119: StokesDOF ***stokes_analytic,***stokes;
1120: StokesDOF stokes_analytic_e[NODES_PER_EL],stokes_e[NODES_PER_EL];
1122: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1123: PetscScalar Ni_p[NODES_PER_EL];
1124: PetscInt ngp;
1125: PetscScalar gp_xi[GAUSS_POINTS][NSD];
1126: PetscScalar gp_weight[GAUSS_POINTS];
1127: PetscInt p,i;
1128: PetscScalar J_p,fac;
1129: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1130: PetscScalar tint_p_ms,tint_p,int_p_ms,int_p;
1131: PetscInt M;
1132: PetscReal xymin[NSD],xymax[NSD];
1135: /* define quadrature rule */
1136: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1138: /* setup for coords */
1139: DMGetCoordinateDM(stokes_da,&cda);
1140: DMGetCoordinatesLocal(stokes_da,&coords);
1141: DMDAVecGetArray(cda,coords,&_coords);
1143: /* setup for analytic */
1144: DMCreateLocalVector(stokes_da,&X_analytic_local);
1145: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1146: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1147: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1149: /* setup for solution */
1150: DMCreateLocalVector(stokes_da,&X_local);
1151: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1152: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1153: DMDAVecGetArray(stokes_da,X_local,&stokes);
1155: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1156: DMGetBoundingBox(stokes_da,xymin,xymax);
1158: h = (xymax[0]-xymin[0])/((PetscReal)(M-1));
1160: tp_L2 = tu_L2 = tu_H1 = 0.0;
1161: tint_p_ms = tint_p = 0.0;
1163: DMDAGetElementsCorners(stokes_da,&sex,&sey,&sez);
1164: DMDAGetElementsSizes(stokes_da,&mx,&my,&mz);
1165: for (ek = sez; ek < sez+mz; ek++) {
1166: for (ej = sey; ej < sey+my; ej++) {
1167: for (ei = sex; ei < sex+mx; ei++) {
1168: /* get coords for the element */
1169: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1170: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1171: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1173: /* evaluate integral */
1174: for (p = 0; p < ngp; p++) {
1175: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1176: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1177: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1178: fac = gp_weight[p]*J_p;
1180: for (i = 0; i < NODES_PER_EL; i++) {
1181: tint_p_ms = tint_p_ms+fac*Ni_p[i]*stokes_analytic_e[i].p_dof;
1182: tint_p = tint_p +fac*Ni_p[i]*stokes_e[i].p_dof;
1183: }
1184: }
1185: }
1186: }
1187: }
1188: MPI_Allreduce(&tint_p_ms,&int_p_ms,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1189: MPI_Allreduce(&tint_p,&int_p,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1191: PetscPrintf(PETSC_COMM_WORLD,"\\int P dv %1.4e (h) %1.4e (ms)\n",(double)PetscRealPart(int_p),(double)PetscRealPart(int_p_ms));
1193: /* remove mine and add manufacture one */
1194: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1195: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1197: {
1198: PetscInt k,L,dof;
1199: PetscScalar *fields;
1200: DMDAGetInfo(stokes_da,0, 0,0,0, 0,0,0, &dof,0,0,0,0,0);
1202: VecGetLocalSize(X_local,&L);
1203: VecGetArray(X_local,&fields);
1204: for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1205: VecRestoreArray(X_local,&fields);
1207: VecGetLocalSize(X,&L);
1208: VecGetArray(X,&fields);
1209: for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1210: VecRestoreArray(X,&fields);
1211: }
1213: DMDAVecGetArray(stokes_da,X_local,&stokes);
1214: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1216: for (ek = sez; ek < sez+mz; ek++) {
1217: for (ej = sey; ej < sey+my; ej++) {
1218: for (ei = sex; ei < sex+mx; ei++) {
1219: /* get coords for the element */
1220: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1221: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1222: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1224: /* evaluate integral */
1225: p_e_L2 = 0.0;
1226: u_e_L2 = 0.0;
1227: u_e_H1 = 0.0;
1228: for (p = 0; p < ngp; p++) {
1229: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1230: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1231: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1232: fac = gp_weight[p]*J_p;
1234: for (i = 0; i < NODES_PER_EL; i++) {
1235: PetscScalar u_error,v_error,w_error;
1237: 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);
1239: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1240: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1241: w_error = stokes_e[i].w_dof-stokes_analytic_e[i].w_dof;
1242: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error+w_error*w_error);
1244: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1245: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error
1246: +GNx_p[2][i]*u_error*GNx_p[2][i]*u_error
1247: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1248: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error
1249: +GNx_p[2][i]*v_error*GNx_p[2][i]*v_error
1250: +GNx_p[0][i]*w_error*GNx_p[0][i]*w_error /* dw/dx */
1251: +GNx_p[1][i]*w_error*GNx_p[1][i]*w_error
1252: +GNx_p[2][i]*w_error*GNx_p[2][i]*w_error);
1253: }
1254: }
1256: tp_L2 += p_e_L2;
1257: tu_L2 += u_e_L2;
1258: tu_H1 += u_e_H1;
1259: }
1260: }
1261: }
1262: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1263: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1264: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1265: p_L2 = PetscSqrtScalar(p_L2);
1266: u_L2 = PetscSqrtScalar(u_L2);
1267: u_H1 = PetscSqrtScalar(u_H1);
1269: 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));
1271: DMDAVecRestoreArray(cda,coords,&_coords);
1273: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1274: VecDestroy(&X_analytic_local);
1275: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1276: VecDestroy(&X_local);
1277: return 0;
1278: }
1280: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da,Vec FIELD,const char file_prefix[])
1281: {
1282: char vtk_filename[PETSC_MAX_PATH_LEN];
1283: PetscMPIInt rank;
1284: MPI_Comm comm;
1285: FILE *vtk_fp = NULL;
1286: const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
1287: PetscInt si,sj,sk,nx,ny,nz,i;
1288: PetscInt f,n_fields,N;
1289: DM cda;
1290: Vec coords;
1291: Vec l_FIELD;
1292: PetscScalar *_L_FIELD;
1293: PetscInt memory_offset;
1294: PetscScalar *buffer;
1298: /* create file name */
1299: PetscObjectGetComm((PetscObject)da,&comm);
1300: MPI_Comm_rank(comm,&rank);
1301: PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"subdomain-%s-p%1.4d.vts",file_prefix,rank);
1303: /* open file and write header */
1304: vtk_fp = fopen(vtk_filename,"w");
1307: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");
1309: /* coords */
1310: DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1311: N = nx * ny * nz;
1313: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n",byte_order);
1314: 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);
1315: 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);
1317: memory_offset = 0;
1319: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <CellData></CellData>\n");
1321: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <Points>\n");
1323: /* copy coordinates */
1324: DMGetCoordinateDM(da,&cda);
1325: DMGetCoordinatesLocal(da,&coords);
1326: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%" PetscInt_FMT "\" />\n",memory_offset);
1327: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N*3;
1329: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </Points>\n");
1331: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PointData Scalars=\" ");
1332: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_fields,0,0,0,0,0);
1333: for (f=0; f<n_fields; f++) {
1334: const char *field_name;
1335: DMDAGetFieldName(da,f,&field_name);
1336: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"%s ",field_name);
1337: }
1338: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\">\n");
1340: for (f=0; f<n_fields; f++) {
1341: const char *field_name;
1343: DMDAGetFieldName(da,f,&field_name);
1344: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <DataArray type=\"Float64\" Name=\"%s\" format=\"appended\" offset=\"%" PetscInt_FMT "\"/>\n", field_name,memory_offset);
1345: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N;
1346: }
1348: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PointData>\n");
1349: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </Piece>\n");
1350: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </StructuredGrid>\n");
1352: PetscMalloc1(N,&buffer);
1353: DMGetLocalVector(da,&l_FIELD);
1354: DMGlobalToLocalBegin(da, FIELD,INSERT_VALUES,l_FIELD);
1355: DMGlobalToLocalEnd(da,FIELD,INSERT_VALUES,l_FIELD);
1356: VecGetArray(l_FIELD,&_L_FIELD);
1358: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <AppendedData encoding=\"raw\">\n");
1359: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"_");
1361: /* write coordinates */
1362: {
1363: int length = sizeof(PetscScalar)*N*3;
1364: PetscScalar *allcoords;
1366: fwrite(&length,sizeof(int),1,vtk_fp);
1367: VecGetArray(coords,&allcoords);
1368: fwrite(allcoords,sizeof(PetscScalar),3*N,vtk_fp);
1369: VecRestoreArray(coords,&allcoords);
1370: }
1371: /* write fields */
1372: for (f=0; f<n_fields; f++) {
1373: int length = sizeof(PetscScalar)*N;
1374: fwrite(&length,sizeof(int),1,vtk_fp);
1375: /* load */
1376: for (i=0; i<N; i++) buffer[i] = _L_FIELD[n_fields*i + f];
1378: /* write */
1379: fwrite(buffer,sizeof(PetscScalar),N,vtk_fp);
1380: }
1381: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\n </AppendedData>\n");
1383: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");
1385: PetscFree(buffer);
1386: VecRestoreArray(l_FIELD,&_L_FIELD);
1387: DMRestoreLocalVector(da,&l_FIELD);
1389: if (vtk_fp) {
1390: fclose(vtk_fp);
1391: vtk_fp = NULL;
1392: }
1394: return 0;
1395: }
1397: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp,PetscInt indent_level,DM da,const char local_file_prefix[])
1398: {
1399: PetscMPIInt size,rank;
1400: MPI_Comm comm;
1401: const PetscInt *lx,*ly,*lz;
1402: PetscInt M,N,P,pM,pN,pP,sum,*olx,*oly,*olz;
1403: PetscInt *osx,*osy,*osz,*oex,*oey,*oez;
1404: PetscInt i,j,k,II,stencil;
1407: /* create file name */
1408: PetscObjectGetComm((PetscObject)da,&comm);
1409: MPI_Comm_size(comm,&size);
1410: MPI_Comm_rank(comm,&rank);
1412: DMDAGetInfo(da,0,&M,&N,&P,&pM,&pN,&pP,0,&stencil,0,0,0,0);
1413: DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
1415: /* generate start,end list */
1416: PetscMalloc1(pM+1,&olx);
1417: PetscMalloc1(pN+1,&oly);
1418: PetscMalloc1(pP+1,&olz);
1419: sum = 0;
1420: for (i=0; i<pM; i++) {
1421: olx[i] = sum;
1422: sum = sum + lx[i];
1423: }
1424: olx[pM] = sum;
1425: sum = 0;
1426: for (i=0; i<pN; i++) {
1427: oly[i] = sum;
1428: sum = sum + ly[i];
1429: }
1430: oly[pN] = sum;
1431: sum = 0;
1432: for (i=0; i<pP; i++) {
1433: olz[i] = sum;
1434: sum = sum + lz[i];
1435: }
1436: olz[pP] = sum;
1438: PetscMalloc1(pM,&osx);
1439: PetscMalloc1(pN,&osy);
1440: PetscMalloc1(pP,&osz);
1441: PetscMalloc1(pM,&oex);
1442: PetscMalloc1(pN,&oey);
1443: PetscMalloc1(pP,&oez);
1444: for (i=0; i<pM; i++) {
1445: osx[i] = olx[i] - stencil;
1446: oex[i] = olx[i] + lx[i] + stencil;
1447: if (osx[i]<0) osx[i]=0;
1448: if (oex[i]>M) oex[i]=M;
1449: }
1451: for (i=0; i<pN; i++) {
1452: osy[i] = oly[i] - stencil;
1453: oey[i] = oly[i] + ly[i] + stencil;
1454: if (osy[i]<0)osy[i]=0;
1455: if (oey[i]>M)oey[i]=N;
1456: }
1457: for (i=0; i<pP; i++) {
1458: osz[i] = olz[i] - stencil;
1459: oez[i] = olz[i] + lz[i] + stencil;
1460: if (osz[i]<0) osz[i]=0;
1461: if (oez[i]>P) oez[i]=P;
1462: }
1464: for (k=0; k<pP; k++) {
1465: for (j=0; j<pN; j++) {
1466: for (i=0; i<pM; i++) {
1467: char name[PETSC_MAX_PATH_LEN];
1468: PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
1469: PetscSNPrintf(name,sizeof(name),"subdomain-%s-p%1.4" PetscInt_FMT ".vts",local_file_prefix,procid);
1470: for (II=0; II<indent_level; II++) PetscFPrintf(PETSC_COMM_SELF,vtk_fp," ");
1472: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<Piece Extent=\"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\" Source=\"%s\"/>\n",
1473: osx[i],oex[i]-1,
1474: osy[j],oey[j]-1,
1475: osz[k],oez[k]-1,name);
1476: }
1477: }
1478: }
1479: PetscFree(olx);
1480: PetscFree(oly);
1481: PetscFree(olz);
1482: PetscFree(osx);
1483: PetscFree(osy);
1484: PetscFree(osz);
1485: PetscFree(oex);
1486: PetscFree(oey);
1487: PetscFree(oez);
1488: return 0;
1489: }
1491: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da,const char file_prefix[],const char local_file_prefix[])
1492: {
1493: MPI_Comm comm;
1494: PetscMPIInt size,rank;
1495: char vtk_filename[PETSC_MAX_PATH_LEN];
1496: FILE *vtk_fp = NULL;
1497: const char *byte_order = PetscBinaryBigEndian() ? "BigEndian" : "LittleEndian";
1498: PetscInt M,N,P,si,sj,sk,nx,ny,nz;
1499: PetscInt i,dofs;
1502: /* only rank-0 generates this file */
1503: PetscObjectGetComm((PetscObject)da,&comm);
1504: MPI_Comm_size(comm,&size);
1505: MPI_Comm_rank(comm,&rank);
1507: if (rank != 0) return 0;
1509: /* create file name */
1510: PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"%s.pvts",file_prefix);
1511: vtk_fp = fopen(vtk_filename,"w");
1514: /* (VTK) generate pvts header */
1515: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");
1517: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"%s\">\n",byte_order);
1519: /* define size of the nodal mesh based on the cell DM */
1520: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,&dofs,0,0,0,0,0);
1521: DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1522: 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 */
1524: /* DUMP THE CELL REFERENCES */
1525: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PCellData>\n");
1526: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PCellData>\n");
1528: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PPoints>\n");
1529: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"%d\"/>\n",NSD);
1530: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PPoints>\n");
1532: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PPointData>\n");
1533: for (i=0; i<dofs; i++) {
1534: const char *fieldname;
1535: DMDAGetFieldName(da,i,&fieldname);
1536: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n",fieldname);
1537: }
1538: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PPointData>\n");
1540: /* write out the parallel information */
1541: DAViewVTK_write_PieceExtend(vtk_fp,2,da,local_file_prefix);
1543: /* close the file */
1544: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PStructuredGrid>\n");
1545: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");
1547: if (vtk_fp) {
1548: fclose(vtk_fp);
1549: vtk_fp = NULL;
1550: }
1551: return 0;
1552: }
1554: PetscErrorCode DAView3DPVTS(DM da, Vec x,const char NAME[])
1555: {
1556: char vts_filename[PETSC_MAX_PATH_LEN];
1557: char pvts_filename[PETSC_MAX_PATH_LEN];
1560: PetscSNPrintf(vts_filename,sizeof(vts_filename),"%s-mesh",NAME);
1561: DAView_3DVTK_StructuredGrid_appended(da,x,vts_filename);
1563: PetscSNPrintf(pvts_filename,sizeof(pvts_filename),"%s-mesh",NAME);
1564: DAView_3DVTK_PStructuredGrid(da,pvts_filename,vts_filename);
1565: return 0;
1566: }
1568: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
1569: {
1570: PetscReal norms[4];
1571: Vec Br,v,w;
1572: Mat A;
1575: KSPGetOperators(ksp,&A,NULL);
1576: MatCreateVecs(A,&w,&v);
1578: KSPBuildResidual(ksp,v,w,&Br);
1580: VecStrideNorm(Br,0,NORM_2,&norms[0]);
1581: VecStrideNorm(Br,1,NORM_2,&norms[1]);
1582: VecStrideNorm(Br,2,NORM_2,&norms[2]);
1583: VecStrideNorm(Br,3,NORM_2,&norms[3]);
1585: VecDestroy(&v);
1586: VecDestroy(&w);
1588: 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]);
1589: return 0;
1590: }
1592: static PetscErrorCode PCMGSetupViaCoarsen(PC pc,DM da_fine)
1593: {
1594: PetscInt nlevels,k;
1595: PETSC_UNUSED PetscInt finest;
1596: DM *da_list,*daclist;
1597: Mat R;
1600: nlevels = 1;
1601: PetscOptionsGetInt(NULL,NULL,"-levels",&nlevels,0);
1603: PetscMalloc1(nlevels,&da_list);
1604: for (k=0; k<nlevels; k++) da_list[k] = NULL;
1605: PetscMalloc1(nlevels,&daclist);
1606: for (k=0; k<nlevels; k++) daclist[k] = NULL;
1608: /* finest grid is nlevels - 1 */
1609: finest = nlevels - 1;
1610: daclist[0] = da_fine;
1611: PetscObjectReference((PetscObject)da_fine);
1612: DMCoarsenHierarchy(da_fine,nlevels-1,&daclist[1]);
1613: for (k=0; k<nlevels; k++) {
1614: da_list[k] = daclist[nlevels-1-k];
1615: DMDASetUniformCoordinates(da_list[k],0.0,1.0,0.0,1.0,0.0,1.0);
1616: }
1618: PCMGSetLevels(pc,nlevels,NULL);
1619: PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
1620: PCMGSetGalerkin(pc,PC_MG_GALERKIN_PMAT);
1622: for (k=1; k<nlevels; k++) {
1623: DMCreateInterpolation(da_list[k-1],da_list[k],&R,NULL);
1624: PCMGSetInterpolation(pc,k,R);
1625: MatDestroy(&R);
1626: }
1628: /* tidy up */
1629: for (k=0; k<nlevels; k++) {
1630: DMDestroy(&da_list[k]);
1631: }
1632: PetscFree(da_list);
1633: PetscFree(daclist);
1634: return 0;
1635: }
1637: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx,PetscInt my,PetscInt mz)
1638: {
1639: DM da_Stokes;
1640: PetscInt u_dof,p_dof,dof,stencil_width;
1641: Mat A,B;
1642: DM vel_cda;
1643: Vec vel_coords;
1644: PetscInt p;
1645: Vec f,X,X1;
1646: DMDACoor3d ***_vel_coords;
1647: PetscInt its;
1648: KSP ksp_S;
1649: PetscInt model_definition = 0;
1650: PetscInt ei,ej,ek,sex,sey,sez,Imx,Jmy,Kmz;
1651: CellProperties cell_properties;
1652: PetscBool write_output = PETSC_FALSE,resolve= PETSC_FALSE;
1655: PetscOptionsGetBool(NULL,NULL,"-resolve",&resolve,NULL);
1656: /* Generate the da for velocity and pressure */
1657: /* Num nodes in each direction is mx+1, my+1, mz+1 */
1658: u_dof = U_DOFS; /* Vx, Vy - velocities */
1659: p_dof = P_DOFS; /* p - pressure */
1660: dof = u_dof+p_dof;
1661: stencil_width = 1;
1662: PetscCall(DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1663: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,NULL,&da_Stokes));
1664: DMSetMatType(da_Stokes,MATAIJ);
1665: DMSetFromOptions(da_Stokes);
1666: DMSetUp(da_Stokes);
1667: DMDASetFieldName(da_Stokes,0,"Vx");
1668: DMDASetFieldName(da_Stokes,1,"Vy");
1669: DMDASetFieldName(da_Stokes,2,"Vz");
1670: DMDASetFieldName(da_Stokes,3,"P");
1672: /* unit box [0,1] x [0,1] x [0,1] */
1673: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.0,1.0);
1675: /* create quadrature point info for PDE coefficients */
1676: CellPropertiesCreate(da_Stokes,&cell_properties);
1678: /* interpolate the coordinates to quadrature points */
1679: DMGetCoordinateDM(da_Stokes,&vel_cda);
1680: DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1681: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1682: DMDAGetElementsCorners(da_Stokes,&sex,&sey,&sez);
1683: DMDAGetElementsSizes(da_Stokes,&Imx,&Jmy,&Kmz);
1684: for (ek = sez; ek < sez+Kmz; ek++) {
1685: for (ej = sey; ej < sey+Jmy; ej++) {
1686: for (ei = sex; ei < sex+Imx; ei++) {
1687: /* get coords for the element */
1688: PetscInt ngp;
1689: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
1690: PetscScalar el_coords[NSD*NODES_PER_EL];
1691: GaussPointCoefficients *cell;
1693: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1694: GetElementCoords3D(_vel_coords,ei,ej,ek,el_coords);
1695: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1697: for (p = 0; p < GAUSS_POINTS; p++) {
1698: PetscScalar xi_p[NSD],Ni_p[NODES_PER_EL];
1699: PetscScalar gp_x,gp_y,gp_z;
1700: PetscInt n;
1702: xi_p[0] = gp_xi[p][0];
1703: xi_p[1] = gp_xi[p][1];
1704: xi_p[2] = gp_xi[p][2];
1705: ShapeFunctionQ13D_Evaluate(xi_p,Ni_p);
1707: gp_x = gp_y = gp_z = 0.0;
1708: for (n = 0; n < NODES_PER_EL; n++) {
1709: gp_x = gp_x+Ni_p[n]*el_coords[NSD*n];
1710: gp_y = gp_y+Ni_p[n]*el_coords[NSD*n+1];
1711: gp_z = gp_z+Ni_p[n]*el_coords[NSD*n+2];
1712: }
1713: cell->gp_coords[NSD*p] = gp_x;
1714: cell->gp_coords[NSD*p+1] = gp_y;
1715: cell->gp_coords[NSD*p+2] = gp_z;
1716: }
1717: }
1718: }
1719: }
1721: PetscOptionsGetInt(NULL,NULL,"-model",&model_definition,NULL);
1723: switch (model_definition) {
1724: case 0: /* isoviscous */
1725: for (ek = sez; ek < sez+Kmz; ek++) {
1726: for (ej = sey; ej < sey+Jmy; ej++) {
1727: for (ei = sex; ei < sex+Imx; ei++) {
1728: GaussPointCoefficients *cell;
1730: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1731: for (p = 0; p < GAUSS_POINTS; p++) {
1732: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1733: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1734: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1736: cell->eta[p] = 1.0;
1738: cell->fx[p] = 0.0*coord_x;
1739: cell->fy[p] = -PetscSinReal(2.2*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1740: cell->fz[p] = 0.0*coord_z;
1741: cell->hc[p] = 0.0;
1742: }
1743: }
1744: }
1745: }
1746: break;
1748: case 1: /* manufactured */
1749: for (ek = sez; ek < sez+Kmz; ek++) {
1750: for (ej = sey; ej < sey+Jmy; ej++) {
1751: for (ei = sex; ei < sex+Imx; ei++) {
1752: PetscReal eta,Fm[NSD],Fc,pos[NSD];
1753: GaussPointCoefficients *cell;
1755: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1756: for (p = 0; p < GAUSS_POINTS; p++) {
1757: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1758: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1759: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1761: pos[0] = coord_x;
1762: pos[1] = coord_y;
1763: pos[2] = coord_z;
1765: evaluate_MS_FrankKamentski(pos,NULL,NULL,&eta,Fm,&Fc);
1766: cell->eta[p] = eta;
1768: cell->fx[p] = Fm[0];
1769: cell->fy[p] = Fm[1];
1770: cell->fz[p] = Fm[2];
1771: cell->hc[p] = Fc;
1772: }
1773: }
1774: }
1775: }
1776: break;
1778: case 2: /* solcx */
1779: for (ek = sez; ek < sez+Kmz; ek++) {
1780: for (ej = sey; ej < sey+Jmy; ej++) {
1781: for (ei = sex; ei < sex+Imx; ei++) {
1782: GaussPointCoefficients *cell;
1784: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1785: for (p = 0; p < GAUSS_POINTS; p++) {
1786: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1787: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1788: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1790: cell->eta[p] = 1.0;
1792: cell->fx[p] = 0.0;
1793: cell->fy[p] = -PetscSinReal(3.0*PETSC_PI*coord_y)*PetscCosReal(1.0*PETSC_PI*coord_x);
1794: cell->fz[p] = 0.0*coord_z;
1795: cell->hc[p] = 0.0;
1796: }
1797: }
1798: }
1799: }
1800: break;
1802: case 3: /* sinker */
1803: for (ek = sez; ek < sez+Kmz; ek++) {
1804: for (ej = sey; ej < sey+Jmy; ej++) {
1805: for (ei = sex; ei < sex+Imx; ei++) {
1806: GaussPointCoefficients *cell;
1808: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1809: for (p = 0; p < GAUSS_POINTS; p++) {
1810: PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p]);
1811: PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1812: PetscReal zp = PetscRealPart(cell->gp_coords[NSD*p+2]);
1814: cell->eta[p] = 1.0e-2;
1815: cell->fx[p] = 0.0;
1816: cell->fy[p] = 0.0;
1817: cell->fz[p] = 0.0;
1818: cell->hc[p] = 0.0;
1820: if ((PetscAbs(xp-0.5) < 0.2) && (PetscAbs(yp-0.5) < 0.2) && (PetscAbs(zp-0.5) < 0.2)) {
1821: cell->eta[p] = 1.0;
1822: cell->fz[p] = 1.0;
1823: }
1825: }
1826: }
1827: }
1828: }
1829: break;
1831: case 4: /* subdomain jumps */
1832: for (ek = sez; ek < sez+Kmz; ek++) {
1833: for (ej = sey; ej < sey+Jmy; ej++) {
1834: for (ei = sex; ei < sex+Imx; ei++) {
1835: PetscReal opts_mag,opts_eta0;
1836: PetscInt px,py,pz;
1837: PetscBool jump;
1838: PetscMPIInt rr;
1839: GaussPointCoefficients *cell;
1841: opts_mag = 1.0;
1842: opts_eta0 = 1.e-2;
1844: PetscOptionsGetReal(NULL,NULL,"-jump_eta0",&opts_eta0,NULL);
1845: PetscOptionsGetReal(NULL,NULL,"-jump_magnitude",&opts_mag,NULL);
1846: DMDAGetInfo(da_Stokes,NULL,NULL,NULL,NULL,&px,&py,&pz,NULL,NULL,NULL,NULL,NULL,NULL);
1847: rr = PetscGlobalRank%(px*py);
1848: if (px%2) jump = (PetscBool)(rr%2);
1849: else jump = (PetscBool)((rr/px)%2 ? rr%2 : !(rr%2));
1850: rr = PetscGlobalRank/(px*py);
1851: if (rr%2) jump = (PetscBool)!jump;
1852: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1853: for (p = 0; p < GAUSS_POINTS; p++) {
1854: PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p]);
1855: PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1857: cell->eta[p] = jump ? PetscPowReal(10.0,opts_mag) : opts_eta0;
1858: cell->fx[p] = 0.0;
1859: cell->fy[p] = -PetscSinReal(2.2*PETSC_PI*yp)*PetscCosReal(1.0*PETSC_PI*xp);
1860: cell->fz[p] = 0.0;
1861: cell->hc[p] = 0.0;
1862: }
1863: }
1864: }
1865: }
1866: break;
1867: default:
1868: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No default model is supported. Choose either -model {0,1,2,3}");
1869: }
1871: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1873: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1874: DMCreateMatrix(da_Stokes,&A);
1875: DMCreateMatrix(da_Stokes,&B);
1876: DMCreateGlobalVector(da_Stokes,&X);
1877: DMCreateGlobalVector(da_Stokes,&f);
1879: /* assemble A11 */
1880: MatZeroEntries(A);
1881: MatZeroEntries(B);
1882: VecZeroEntries(f);
1884: AssembleA_Stokes(A,da_Stokes,cell_properties);
1885: MatViewFromOptions(A,NULL,"-amat_view");
1886: AssembleA_PCStokes(B,da_Stokes,cell_properties);
1887: MatViewFromOptions(B,NULL,"-bmat_view");
1888: /* build force vector */
1889: AssembleF_Stokes(f,da_Stokes,cell_properties);
1891: /* SOLVE */
1892: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1893: KSPSetOptionsPrefix(ksp_S,"stokes_");
1894: KSPSetOperators(ksp_S,A,B);
1895: KSPSetFromOptions(ksp_S);
1897: {
1898: PC pc;
1899: const PetscInt ufields[] = {0,1,2},pfields[] = {3};
1900: KSPGetPC(ksp_S,&pc);
1901: PCFieldSplitSetBlockSize(pc,4);
1902: PCFieldSplitSetFields(pc,"u",3,ufields,ufields);
1903: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1904: }
1906: {
1907: PC pc;
1908: PetscBool same = PETSC_FALSE;
1909: KSPGetPC(ksp_S,&pc);
1910: PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
1911: if (same) PCMGSetupViaCoarsen(pc,da_Stokes);
1912: }
1914: {
1915: PC pc;
1916: PetscBool same = PETSC_FALSE;
1917: KSPGetPC(ksp_S,&pc);
1918: PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&same);
1919: if (same) KSPSetOperators(ksp_S,A,A);
1920: }
1922: {
1923: PetscBool stokes_monitor = PETSC_FALSE;
1924: PetscOptionsGetBool(NULL,NULL,"-stokes_ksp_monitor_blocks",&stokes_monitor,0);
1925: if (stokes_monitor) KSPMonitorSet(ksp_S,KSPMonitorStokesBlocks,NULL,NULL);
1926: }
1928: if (resolve) {
1929: /* Test changing matrix data structure and resolve */
1930: VecDuplicate(X,&X1);
1931: }
1933: KSPSolve(ksp_S,f,X);
1934: if (resolve) {
1935: Mat C;
1936: MatDuplicate(A,MAT_COPY_VALUES,&C);
1937: KSPSetOperators(ksp_S,C,C);
1938: KSPSolve(ksp_S,f,X1);
1939: MatDestroy(&C);
1940: VecDestroy(&X1);
1941: }
1943: PetscOptionsGetBool(NULL,NULL,"-write_pvts",&write_output,NULL);
1944: if (write_output) DAView3DPVTS(da_Stokes,X,"up");
1945: {
1946: PetscBool flg = PETSC_FALSE;
1947: char filename[PETSC_MAX_PATH_LEN];
1948: PetscOptionsGetString(NULL,NULL,"-write_binary",filename,sizeof(filename),&flg);
1949: if (flg) {
1950: PetscViewer viewer;
1951: /* PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer); */
1952: PetscViewerVTKOpen(PETSC_COMM_WORLD,"ex42.vts",FILE_MODE_WRITE,&viewer);
1953: VecView(X,viewer);
1954: PetscViewerDestroy(&viewer);
1955: }
1956: }
1957: KSPGetIterationNumber(ksp_S,&its);
1959: /* verify */
1960: if (model_definition == 1) {
1961: DM da_Stokes_analytic;
1962: Vec X_analytic;
1964: DMDACreateManufacturedSolution(mx,my,mz,&da_Stokes_analytic,&X_analytic);
1965: if (write_output) {
1966: DAView3DPVTS(da_Stokes_analytic,X_analytic,"ms");
1967: }
1968: DMDAIntegrateErrors3D(da_Stokes_analytic,X,X_analytic);
1969: if (write_output) {
1970: DAView3DPVTS(da_Stokes,X,"up2");
1971: }
1972: DMDestroy(&da_Stokes_analytic);
1973: VecDestroy(&X_analytic);
1974: }
1976: KSPDestroy(&ksp_S);
1977: VecDestroy(&X);
1978: VecDestroy(&f);
1979: MatDestroy(&A);
1980: MatDestroy(&B);
1982: CellPropertiesDestroy(&cell_properties);
1983: DMDestroy(&da_Stokes);
1984: return 0;
1985: }
1987: int main(int argc,char **args)
1988: {
1989: PetscInt mx,my,mz;
1992: PetscInitialize(&argc,&args,(char*)0,help);
1993: mx = my = mz = 10;
1994: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1995: my = mx; mz = mx;
1996: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1997: PetscOptionsGetInt(NULL,NULL,"-mz",&mz,NULL);
1998: solve_stokes_3d_coupled(mx,my,mz);
1999: PetscFinalize();
2000: return 0;
2001: }
2003: /*TEST
2005: test:
2006: args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type lu
2008: test:
2009: suffix: 2
2010: nsize: 3
2011: args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type redundant -stokes_redundant_pc_type lu
2013: test:
2014: suffix: bddc_stokes
2015: nsize: 6
2016: 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
2018: test:
2019: suffix: bddc_stokes_deluxe
2020: nsize: 6
2021: 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
2023: test:
2024: requires: !single
2025: suffix: bddc_stokes_subdomainjump_deluxe
2026: nsize: 9
2027: 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
2029: test:
2030: requires: !single
2031: suffix: 3
2032: args: -stokes_ksp_converged_reason -stokes_pc_type fieldsplit -resolve
2034: test:
2035: suffix: tut_1
2036: nsize: 4
2037: requires: !single
2038: args: -stokes_ksp_monitor
2040: test:
2041: suffix: tut_2
2042: nsize: 4
2043: requires: !single
2044: args: -stokes_ksp_monitor -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type schur
2046: test:
2047: suffix: tut_3
2048: nsize: 4
2049: requires: !single
2050: args: -mx 20 -stokes_ksp_monitor -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_type schur
2052: TEST*/