Actual source code: ex42.c
petsc-3.3-p7 2013-05-11
1: static char help[] =
2: "Solves the incompressible, variable viscosity stokes equation in 3d using Q1Q1 elements, \n\
3: stabilized with Bochev's polynomial projection method. Note that implementation here assumes \n\
4: all boundaries are free-slip, i.e. zero normal flow and zero tangential stress \n\
5: -mx : number elements in x-direction \n\
6: -my : number elements in y-direction \n\
7: -mz : number elements in z-direction \n\
8: -stokes_ksp_monitor_blocks : active monitor for each component u,v,w,p \n\
9: -model : defines viscosity and forcing function. Choose either: 0 (isoviscous), 1 (manufactured-broken), 2 (solcx-2d), 3 (sinker) \n";
11: /* Contributed by Dave May */
13: #include petscksp.h
14: #include petscpcmg.h
15: #include petscdmda.h
17: #define PROFILE_TIMING
18: #define ASSEMBLE_LOWER_TRIANGULAR
20: #define NSD 3 /* number of spatial dimensions */
21: #define NODES_PER_EL 8 /* nodes per element */
22: #define U_DOFS 3 /* degrees of freedom per velocity node */
23: #define P_DOFS 1 /* degrees of freedom per pressure node */
24: #define GAUSS_POINTS 8
26: /* Gauss point based evaluation */
27: typedef struct {
28: PetscScalar gp_coords[NSD*GAUSS_POINTS];
29: PetscScalar eta[GAUSS_POINTS];
30: PetscScalar fx[GAUSS_POINTS];
31: PetscScalar fy[GAUSS_POINTS];
32: PetscScalar fz[GAUSS_POINTS];
33: PetscScalar hc[GAUSS_POINTS];
34: } GaussPointCoefficients;
36: typedef struct {
37: PetscScalar u_dof;
38: PetscScalar v_dof;
39: PetscScalar w_dof;
40: PetscScalar p_dof;
41: } StokesDOF;
43: typedef struct _p_CellProperties *CellProperties;
44: struct _p_CellProperties {
45: PetscInt ncells;
46: PetscInt mx,my,mz;
47: PetscInt sex,sey,sez;
48: GaussPointCoefficients *gpc;
49: };
51: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz);
53: /* elements */
56: PetscErrorCode CellPropertiesCreate(DM da_stokes,CellProperties *C)
57: {
59: CellProperties cells;
60: PetscInt mx,my,mz,sex,sey,sez;
61:
63: PetscMalloc(sizeof(struct _p_CellProperties),&cells);
64:
65: DMDAGetElementCorners(da_stokes,&sex,&sey,&sez,&mx,&my,&mz);
66: cells->mx = mx;
67: cells->my = my;
68: cells->mz = mz;
69: cells->ncells = mx * my * mz;
70: cells->sex = sex;
71: cells->sey = sey;
72: cells->sez = sez;
73: PetscMalloc(sizeof(GaussPointCoefficients)*mx*my*mz,&cells->gpc);
74:
75: *C = cells;
76: return(0);
77: }
81: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
82: {
84: CellProperties cells;
87: if (!C) { return(0); }
88: cells = *C;
89: PetscFree(cells->gpc);
90: PetscFree(cells);
91: *C = PETSC_NULL;
92: return(0);
93: }
97: PetscErrorCode CellPropertiesGetCell(CellProperties C,PetscInt I,PetscInt J,PetscInt K,GaussPointCoefficients **G)
98: {
100: *G = &C->gpc[ (I-C->sex) + (J-C->sey)*C->mx + (K-C->sez)*C->mx*C->my ];
101: return(0);
102: }
104: /* FEM routines */
105: /*
106: Element: Local basis function ordering
107: 1-----2
108: | |
109: | |
110: 0-----3
111: */
112: static void ShapeFunctionQ13D_Evaluate(PetscScalar _xi[],PetscScalar Ni[])
113: {
114: PetscReal xi = PetscRealPart(_xi[0]);
115: PetscReal eta = PetscRealPart(_xi[1]);
116: PetscReal zeta = PetscRealPart(_xi[2]);
118: Ni[0] = 0.125 * ( 1.0 - xi ) * ( 1.0 - eta ) * ( 1.0 - zeta );
119: Ni[1] = 0.125 * ( 1.0 - xi ) * ( 1.0 + eta ) * ( 1.0 - zeta );
120: Ni[2] = 0.125 * ( 1.0 + xi ) * ( 1.0 + eta ) * ( 1.0 - zeta );
121: Ni[3] = 0.125 * ( 1.0 + xi ) * ( 1.0 - eta ) * ( 1.0 - zeta );
123: Ni[4] = 0.125 * ( 1.0 - xi ) * ( 1.0 - eta ) * ( 1.0 + zeta );
124: Ni[5] = 0.125 * ( 1.0 - xi ) * ( 1.0 + eta ) * ( 1.0 + zeta );
125: Ni[6] = 0.125 * ( 1.0 + xi ) * ( 1.0 + eta ) * ( 1.0 + zeta );
126: Ni[7] = 0.125 * ( 1.0 + xi ) * ( 1.0 - eta ) * ( 1.0 + zeta );
127: }
129: static void ShapeFunctionQ13D_Evaluate_dxi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
130: {
131: PetscReal xi = PetscRealPart(_xi[0]);
132: PetscReal eta = PetscRealPart(_xi[1]);
133: PetscReal zeta = PetscRealPart(_xi[2]);
134: /* xi deriv */
135: GNi[0][0] = - 0.125 * ( 1.0 - eta ) * ( 1.0 - zeta );
136: GNi[0][1] = - 0.125 * ( 1.0 + eta ) * ( 1.0 - zeta );
137: GNi[0][2] = 0.125 * ( 1.0 + eta ) * ( 1.0 - zeta );
138: GNi[0][3] = 0.125 * ( 1.0 - eta ) * ( 1.0 - zeta );
140: GNi[0][4] = - 0.125 * ( 1.0 - eta ) * ( 1.0 + zeta );
141: GNi[0][5] = - 0.125 * ( 1.0 + eta ) * ( 1.0 + zeta );
142: GNi[0][6] = 0.125 * ( 1.0 + eta ) * ( 1.0 + zeta );
143: GNi[0][7] = 0.125 * ( 1.0 - eta ) * ( 1.0 + zeta );
144: /* eta deriv */
145: GNi[1][0] = - 0.125 * ( 1.0 - xi ) * ( 1.0 - zeta );
146: GNi[1][1] = 0.125 * ( 1.0 - xi ) * ( 1.0 - zeta );
147: GNi[1][2] = 0.125 * ( 1.0 + xi ) * ( 1.0 - zeta );
148: GNi[1][3] = - 0.125 * ( 1.0 + xi ) * ( 1.0 - zeta );
150: GNi[1][4] = - 0.125 * ( 1.0 - xi ) * ( 1.0 + zeta );
151: GNi[1][5] = 0.125 * ( 1.0 - xi ) * ( 1.0 + zeta );
152: GNi[1][6] = 0.125 * ( 1.0 + xi ) * ( 1.0 + zeta );
153: GNi[1][7] = - 0.125 * ( 1.0 + xi ) * ( 1.0 + zeta );
154: /* zeta deriv */
155: GNi[2][0] = -0.125 * ( 1.0 - xi ) * ( 1.0 - eta );
156: GNi[2][1] = -0.125 * ( 1.0 - xi ) * ( 1.0 + eta );
157: GNi[2][2] = -0.125 * ( 1.0 + xi ) * ( 1.0 + eta );
158: GNi[2][3] = -0.125 * ( 1.0 + xi ) * ( 1.0 - eta );
160: GNi[2][4] = 0.125 * ( 1.0 - xi ) * ( 1.0 - eta );
161: GNi[2][5] = 0.125 * ( 1.0 - xi ) * ( 1.0 + eta );
162: GNi[2][6] = 0.125 * ( 1.0 + xi ) * ( 1.0 + eta );
163: GNi[2][7] = 0.125 * ( 1.0 + xi ) * ( 1.0 - eta );
164: }
166: static void matrix_inverse_3x3(PetscScalar A[3][3],PetscScalar B[3][3])
167: {
168: PetscScalar t4, t6, t8, t10, t12, t14, t17;
170: t4 = A[2][0] * A[0][1];
171: t6 = A[2][0] * A[0][2];
172: t8 = A[1][0] * A[0][1];
173: t10 = A[1][0] * A[0][2];
174: t12 = A[0][0] * A[1][1];
175: t14 = A[0][0] * A[1][2];
176: 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]);
178: B[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) * t17;
179: B[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) * t17;
180: B[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) * t17;
181: B[1][0] = -(-A[2][0] * A[1][2] + A[1][0] * A[2][2]) * t17;
182: B[1][1] = (-t6 + A[0][0] * A[2][2]) * t17;
183: B[1][2] = -(-t10 + t14) * t17;
184: B[2][0] = (-A[2][0] * A[1][1] + A[1][0] * A[2][1]) * t17;
185: B[2][1] = -(-t4 + A[0][0] * A[2][1]) * t17;
186: B[2][2] = (-t8 + t12) * t17;
187: }
189: static void ShapeFunctionQ13D_Evaluate_dx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
190: {
191: PetscScalar J00,J01,J02,J10,J11,J12,J20,J21,J22;
192: PetscInt n;
193: PetscScalar iJ[3][3],JJ[3][3];
195: J00 = J01 = J02 = 0.0;
196: J10 = J11 = J12 = 0.0;
197: J20 = J21 = J22 = 0.0;
198: for (n=0; n<NODES_PER_EL; n++) {
199: PetscScalar cx = coords[ NSD*n + 0 ];
200: PetscScalar cy = coords[ NSD*n + 1 ];
201: PetscScalar cz = coords[ NSD*n + 2 ];
203: /* J_ij = d(x_j) / d(xi_i) *//* J_ij = \sum _I GNi[j][I} * x_i */
204: J00 = J00 + GNi[0][n] * cx; /* J_xx */
205: J01 = J01 + GNi[0][n] * cy; /* J_xy = dx/deta */
206: J02 = J02 + GNi[0][n] * cz; /* J_xz = dx/dzeta */
208: J10 = J10 + GNi[1][n] * cx; /* J_yx = dy/dxi */
209: J11 = J11 + GNi[1][n] * cy; /* J_yy */
210: J12 = J12 + GNi[1][n] * cz; /* J_yz */
212: J20 = J20 + GNi[2][n] * cx; /* J_zx */
213: J21 = J21 + GNi[2][n] * cy; /* J_zy */
214: J22 = J22 + GNi[2][n] * cz; /* J_zz */
215: }
217: JJ[0][0] = J00; JJ[0][1] = J01; JJ[0][2] = J02;
218: JJ[1][0] = J10; JJ[1][1] = J11; JJ[1][2] = J12;
219: JJ[2][0] = J20; JJ[2][1] = J21; JJ[2][2] = J22;
221: matrix_inverse_3x3(JJ,iJ);
223: *det_J = J00*J11*J22 - J00*J12*J21 - J10*J01*J22 + J10*J02*J21 + J20*J01*J12 - J20*J02*J11;
225: for (n=0; n<NODES_PER_EL; n++) {
226: GNx[0][n] = GNi[0][n]*iJ[0][0] + GNi[1][n]*iJ[0][1] + GNi[2][n]*iJ[0][2];
227: GNx[1][n] = GNi[0][n]*iJ[1][0] + GNi[1][n]*iJ[1][1] + GNi[2][n]*iJ[1][2];
228: GNx[2][n] = GNi[0][n]*iJ[2][0] + GNi[1][n]*iJ[2][1] + GNi[2][n]*iJ[2][2];
229: }
230: }
232: static void ConstructGaussQuadrature3D(PetscInt *ngp,PetscScalar gp_xi[][NSD],PetscScalar gp_weight[])
233: {
234: *ngp = 8;
235: gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919; gp_xi[0][2] = -0.57735026919;
236: gp_xi[1][0] = -0.57735026919; gp_xi[1][1] = 0.57735026919; gp_xi[1][2] = -0.57735026919;
237: gp_xi[2][0] = 0.57735026919; gp_xi[2][1] = 0.57735026919; gp_xi[2][2] = -0.57735026919;
238: gp_xi[3][0] = 0.57735026919; gp_xi[3][1] = -0.57735026919; gp_xi[3][2] = -0.57735026919;
240: gp_xi[4][0] = -0.57735026919; gp_xi[4][1] = -0.57735026919; gp_xi[4][2] = 0.57735026919;
241: gp_xi[5][0] = -0.57735026919; gp_xi[5][1] = 0.57735026919; gp_xi[5][2] = 0.57735026919;
242: gp_xi[6][0] = 0.57735026919; gp_xi[6][1] = 0.57735026919; gp_xi[6][2] = 0.57735026919;
243: gp_xi[7][0] = 0.57735026919; gp_xi[7][1] = -0.57735026919; gp_xi[7][2] = 0.57735026919;
245: gp_weight[0] = 1.0;
246: gp_weight[1] = 1.0;
247: gp_weight[2] = 1.0;
248: gp_weight[3] = 1.0;
250: gp_weight[4] = 1.0;
251: gp_weight[5] = 1.0;
252: gp_weight[6] = 1.0;
253: gp_weight[7] = 1.0;
254: }
256: /* procs to the left claim the ghost node as their element */
259: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
260: {
261: PetscInt m,n,p,M,N,P;
262: PetscInt sx,sy,sz;
266: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
267: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
269: if (mxl) {
270: *mxl = m;
271: if ((sx+m) == M) { /* last proc */
272: *mxl = m-1;
273: }
274: }
275: if (myl) {
276: *myl = n;
277: if ((sy+n) == N) { /* last proc */
278: *myl = n-1;
279: }
280: }
281: if (mzl) {
282: *mzl = p;
283: if ((sz+p) == P) { /* last proc */
284: *mzl = p-1;
285: }
286: }
287: return(0);
288: }
292: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
293: {
294: PetscInt si,sj,sk;
298: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
300: if (sx) {
301: *sx = si;
302: if (si != 0) {
303: *sx = si+1;
304: }
305: }
306: if (sy) {
307: *sy = sj;
308: if (sj != 0) {
309: *sy = sj+1;
310: }
311: }
312: if (sz) {
313: *sz = sk;
314: if (sk != 0) {
315: *sz = sk+1;
316: }
317: }
318: DMDAGetLocalElementSize(da,mx,my,mz);
319: return(0);
320: }
322: /*
323: i,j are the element indices
324: The unknown is a vector quantity.
325: The s[].c is used to indicate the degree of freedom.
326: */
329: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j,PetscInt k)
330: {
331: PetscInt n;
333: /* velocity */
334: n = 0;
335: /* node 0 */
336: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++; /* Vx0 */
337: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++; /* Vy0 */
338: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++; /* Vz0 */
340: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
341: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
342: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
344: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
345: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
346: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
348: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++;
349: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++;
350: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++;
352: /* */
353: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++; /* Vx4 */
354: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++; /* Vy4 */
355: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++; /* Vz4 */
357: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
358: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
359: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
361: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
362: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
363: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
365: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++;
366: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++;
367: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++;
369: /* pressure */
370: n = 0;
371: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++; /* P0 */
372: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
373: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
374: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++;
376: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++; /* P0 */
377: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
378: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
379: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++;
380: return(0);
381: }
385: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords,PetscInt i,PetscInt j,PetscInt k,PetscScalar el_coord[])
386: {
388: /* get coords for the element */
389: el_coord[0] = coords[k ][j ][i ].x;
390: el_coord[1] = coords[k ][j ][i ].y;
391: el_coord[2] = coords[k ][j ][i ].z;
393: el_coord[3] = coords[k ][j+1][i].x;
394: el_coord[4] = coords[k ][j+1][i].y;
395: el_coord[5] = coords[k ][j+1][i].z;
397: el_coord[6] = coords[k ][j+1][i+1].x;
398: el_coord[7] = coords[k ][j+1][i+1].y;
399: el_coord[8] = coords[k ][j+1][i+1].z;
401: el_coord[9 ] = coords[k ][j ][i+1].x;
402: el_coord[10] = coords[k ][j ][i+1].y;
403: el_coord[11] = coords[k ][j ][i+1].z;
405: el_coord[12] = coords[k+1][j ][i ].x;
406: el_coord[13] = coords[k+1][j ][i ].y;
407: el_coord[14] = coords[k+1][j ][i ].z;
409: el_coord[15] = coords[k+1][j+1][i ].x;
410: el_coord[16] = coords[k+1][j+1][i ].y;
411: el_coord[17] = coords[k+1][j+1][i ].z;
413: el_coord[18] = coords[k+1][j+1][i+1].x;
414: el_coord[19] = coords[k+1][j+1][i+1].y;
415: el_coord[20] = coords[k+1][j+1][i+1].z;
417: el_coord[21] = coords[k+1][j ][i+1].x;
418: el_coord[22] = coords[k+1][j ][i+1].y;
419: el_coord[23] = coords[k+1][j ][i+1].z;
420: return(0);
421: }
425: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field,PetscInt i,PetscInt j,PetscInt k,StokesDOF nodal_fields[])
426: {
428: /* get the nodal fields for u */
429: nodal_fields[0].u_dof = field[k ][j ][i ].u_dof;
430: nodal_fields[0].v_dof = field[k ][j ][i ].v_dof;
431: nodal_fields[0].w_dof = field[k ][j ][i ].w_dof;
433: nodal_fields[1].u_dof = field[k ][j+1][i ].u_dof;
434: nodal_fields[1].v_dof = field[k ][j+1][i ].v_dof;
435: nodal_fields[1].w_dof = field[k ][j+1][i ].w_dof;
437: nodal_fields[2].u_dof = field[k ][j+1][i+1].u_dof;
438: nodal_fields[2].v_dof = field[k ][j+1][i+1].v_dof;
439: nodal_fields[2].w_dof = field[k ][j+1][i+1].w_dof;
441: nodal_fields[3].u_dof = field[k ][j ][i+1].u_dof;
442: nodal_fields[3].v_dof = field[k ][j ][i+1].v_dof;
443: nodal_fields[3].w_dof = field[k ][j ][i+1].w_dof;
445: nodal_fields[4].u_dof = field[k+1][j ][i ].u_dof;
446: nodal_fields[4].v_dof = field[k+1][j ][i ].v_dof;
447: nodal_fields[4].w_dof = field[k+1][j ][i ].w_dof;
449: nodal_fields[5].u_dof = field[k+1][j+1][i ].u_dof;
450: nodal_fields[5].v_dof = field[k+1][j+1][i ].v_dof;
451: nodal_fields[5].w_dof = field[k+1][j+1][i ].w_dof;
453: nodal_fields[6].u_dof = field[k+1][j+1][i+1].u_dof;
454: nodal_fields[6].v_dof = field[k+1][j+1][i+1].v_dof;
455: nodal_fields[6].w_dof = field[k+1][j+1][i+1].w_dof;
457: nodal_fields[7].u_dof = field[k+1][j ][i+1].u_dof;
458: nodal_fields[7].v_dof = field[k+1][j ][i+1].v_dof;
459: nodal_fields[7].w_dof = field[k+1][j ][i+1].w_dof;
461: /* get the nodal fields for p */
462: nodal_fields[0].p_dof = field[k ][j ][i ].p_dof;
463: nodal_fields[1].p_dof = field[k ][j+1][i ].p_dof;
464: nodal_fields[2].p_dof = field[k ][j+1][i+1].p_dof;
465: nodal_fields[3].p_dof = field[k ][j ][i+1].p_dof;
467: nodal_fields[4].p_dof = field[k+1][j ][i ].p_dof;
468: nodal_fields[5].p_dof = field[k+1][j+1][i ].p_dof;
469: nodal_fields[6].p_dof = field[k+1][j+1][i+1].p_dof;
470: nodal_fields[7].p_dof = field[k+1][j ][i+1].p_dof;
471: return(0);
472: }
474: static PetscInt ASS_MAP_wIwDI_uJuDJ(
475: PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,
476: PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)
477: {
478: PetscInt ij;
479: PETSC_UNUSED PetscInt r,c,nr,nc;
481: nr = w_NPE*w_dof;
482: nc = u_NPE*u_dof;
484: r = w_dof*wi+wd;
485: c = u_dof*ui+ud;
487: ij = r*nc+c;
489: return ij;
490: }
494: static PetscErrorCode DMDASetValuesLocalStencil3D_ADD_VALUES(StokesDOF ***fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
495: {
496: PetscInt n,I,J,K;
499: for (n = 0; n<NODES_PER_EL; n++) {
500: I = u_eqn[NSD*n ].i;
501: J = u_eqn[NSD*n ].j;
502: K = u_eqn[NSD*n ].k;
503: fields_F[K][J][I].u_dof = fields_F[K][J][I].u_dof+Fe_u[NSD*n ];
505: I = u_eqn[NSD*n+1].i;
506: J = u_eqn[NSD*n+1].j;
507: K = u_eqn[NSD*n+1].k;
508: fields_F[K][J][I].v_dof = fields_F[K][J][I].v_dof+Fe_u[NSD*n+1];
510: I = u_eqn[NSD*n+2].i;
511: J = u_eqn[NSD*n+2].j;
512: K = u_eqn[NSD*n+2].k;
513: fields_F[K][J][I].w_dof = fields_F[K][J][I].w_dof+Fe_u[NSD*n+2];
515: I = p_eqn[n].i;
516: J = p_eqn[n].j;
517: K = p_eqn[n].k;
518: fields_F[K][J][I].p_dof = fields_F[K][J][I].p_dof+Fe_p[n ];
520: }
521: return(0);
522: }
524: static void FormStressOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
525: {
526: PetscInt ngp;
527: PetscScalar gp_xi[GAUSS_POINTS][NSD];
528: PetscScalar gp_weight[GAUSS_POINTS];
529: PetscInt p,i,j,k;
530: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
531: PetscScalar J_p,tildeD[6];
532: PetscScalar B[6][U_DOFS*NODES_PER_EL];
533: const PetscInt nvdof = U_DOFS*NODES_PER_EL;
535: /* define quadrature rule */
536: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
538: /* evaluate integral */
539: for (p = 0; p < ngp; p++) {
540: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
541: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
543: for (i = 0; i < NODES_PER_EL; i++) {
544: PetscScalar d_dx_i = GNx_p[0][i];
545: PetscScalar d_dy_i = GNx_p[1][i];
546: PetscScalar d_dz_i = GNx_p[2][i];
548: B[0][3*i ] = d_dx_i; B[0][3*i+1] = 0.0; B[0][3*i+2] = 0.0;
549: B[1][3*i ] = 0.0; B[1][3*i+1] = d_dy_i; B[1][3*i+2] = 0.0;
550: B[2][3*i ] = 0.0; B[2][3*i+1] = 0.0; B[2][3*i+2] = d_dz_i;
552: B[3][3*i] = d_dy_i; B[3][3*i+1] = d_dx_i; B[3][3*i+2] = 0.0; /* e_xy */
553: B[4][3*i] = d_dz_i; B[4][3*i+1] = 0.0; B[4][3*i+2] = d_dx_i;/* e_xz */
554: B[5][3*i] = 0.0; B[5][3*i+1] = d_dz_i; B[5][3*i+2] = d_dy_i;/* e_yz */
555: }
558: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
559: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
560: tildeD[2] = 2.0*gp_weight[p]*J_p*eta[p];
562: tildeD[3] = gp_weight[p]*J_p*eta[p];
563: tildeD[4] = gp_weight[p]*J_p*eta[p];
564: tildeD[5] = gp_weight[p]*J_p*eta[p];
566: /* form Bt tildeD B */
567: /*
568: Ke_ij = Bt_ik . D_kl . B_lj
569: = B_ki . D_kl . B_lj
570: = B_ki . D_kk . B_kj
571: */
572: for (i = 0; i < nvdof; i++) {
573: for (j = i; j < nvdof; j++) {
574: for (k = 0; k < 6; k++) {
575: Ke[i*nvdof+j] += B[k][i]*tildeD[k]*B[k][j];
576: }
577: }
578: }
580: }
581: /* fill lower triangular part */
582: #ifdef ASSEMBLE_LOWER_TRIANGULAR
583: for (i = 0; i < nvdof; i++) {
584: for (j = i; j < nvdof; j++) {
585: Ke[j*nvdof+i] = Ke[i*nvdof+j];
586: }
587: }
588: #endif
589: }
591: static void FormGradientOperatorQ13D(PetscScalar Ke[],PetscScalar coords[])
592: {
593: PetscInt ngp;
594: PetscScalar gp_xi[GAUSS_POINTS][NSD];
595: PetscScalar gp_weight[GAUSS_POINTS];
596: PetscInt p,i,j,di;
597: PetscScalar Ni_p[NODES_PER_EL];
598: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
599: PetscScalar J_p,fac;
601: /* define quadrature rule */
602: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
604: /* evaluate integral */
605: for (p = 0; p < ngp; p++) {
606: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
607: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
608: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
609: fac = gp_weight[p]*J_p;
611: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
612: for (di = 0; di < NSD; di++) { /* u dofs */
613: for (j = 0; j < NODES_PER_EL; j++) { /* p nodes, p dofs = 1 (ie no loop) */
614: PetscInt IJ;
615: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,3,j,0,NODES_PER_EL,1);
617: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
618: }
619: }
620: }
621: }
622: }
624: static void FormDivergenceOperatorQ13D(PetscScalar De[],PetscScalar coords[])
625: {
626: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
627: PetscInt i,j;
628: PetscInt nr_g,nc_g;
630: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
631: FormGradientOperatorQ13D(Ge,coords);
633: nr_g = U_DOFS*NODES_PER_EL;
634: nc_g = P_DOFS*NODES_PER_EL;
636: for (i = 0; i < nr_g; i++) {
637: for (j = 0; j < nc_g; j++) {
638: De[nr_g*j+i] = Ge[nc_g*i+j];
639: }
640: }
641: }
643: static void FormStabilisationOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
644: {
645: PetscInt ngp;
646: PetscScalar gp_xi[GAUSS_POINTS][NSD];
647: PetscScalar gp_weight[GAUSS_POINTS];
648: PetscInt p,i,j;
649: PetscScalar Ni_p[NODES_PER_EL];
650: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
651: PetscScalar J_p,fac,eta_avg;
653: /* define quadrature rule */
654: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
656: /* evaluate integral */
657: for (p = 0; p < ngp; p++) {
658: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
659: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
660: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
661: fac = gp_weight[p]*J_p;
662: /*
663: for (i = 0; i < NODES_PER_EL; i++) {
664: for (j = i; j < NODES_PER_EL; j++) {
665: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
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*(Ni_p[i]*Ni_p[j]-0.015625);
673: }
674: }
675: }
677: /* scale */
678: eta_avg = 0.0;
679: for (p = 0; p < ngp; p++) {
680: eta_avg += eta[p];
681: }
682: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
683: fac = 1.0/eta_avg;
685: /*
686: for (i = 0; i < NODES_PER_EL; i++) {
687: for (j = i; j < NODES_PER_EL; j++) {
688: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
689: #ifdef ASSEMBLE_LOWER_TRIANGULAR
690: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
691: #endif
692: }
693: }
694: */
696: for (i = 0; i < NODES_PER_EL; i++) {
697: for (j = 0; j < NODES_PER_EL; j++) {
698: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
699: }
700: }
701: }
703: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
704: {
705: PetscInt ngp;
706: PetscScalar gp_xi[GAUSS_POINTS][NSD];
707: PetscScalar gp_weight[GAUSS_POINTS];
708: PetscInt p,i,j;
709: PetscScalar Ni_p[NODES_PER_EL];
710: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
711: PetscScalar J_p,fac,eta_avg;
713: /* define quadrature rule */
714: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
716: /* evaluate integral */
717: for (p = 0; p < ngp; p++) {
718: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
719: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
720: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
721: fac = gp_weight[p]*J_p;
723: /*
724: for (i = 0; i < NODES_PER_EL; i++) {
725: for (j = i; j < NODES_PER_EL; j++) {
726: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
727: }
728: }
729: */
731: for (i = 0; i < NODES_PER_EL; i++) {
732: for (j = 0; j < NODES_PER_EL; j++) {
733: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
734: }
735: }
736: }
738: /* scale */
739: eta_avg = 0.0;
740: for (p = 0; p < ngp; p++) {
741: eta_avg += eta[p];
742: }
743: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
744: fac = 1.0/eta_avg;
745: /*
746: for (i = 0; i < NODES_PER_EL; i++) {
747: for (j = i; j < NODES_PER_EL; j++) {
748: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
749: #ifdef ASSEMBLE_LOWER_TRIANGULAR
750: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
751: #endif
752: }
753: }
754: */
756: for (i = 0; i < NODES_PER_EL; i++) {
757: for (j = 0; j < NODES_PER_EL; j++) {
758: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
759: }
760: }
761: }
763: static void FormMomentumRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[],PetscScalar fz[])
764: {
765: PetscInt ngp;
766: PetscScalar gp_xi[GAUSS_POINTS][NSD];
767: PetscScalar gp_weight[GAUSS_POINTS];
768: PetscInt p,i;
769: PetscScalar Ni_p[NODES_PER_EL];
770: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
771: PetscScalar J_p,fac;
773: /* define quadrature rule */
774: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
776: /* evaluate integral */
777: for (p = 0; p < ngp; p++) {
778: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
779: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
780: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
781: fac = gp_weight[p]*J_p;
783: for (i = 0; i < NODES_PER_EL; i++) {
784: Fe[NSD*i ] -= fac*Ni_p[i]*fx[p];
785: Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
786: Fe[NSD*i+2] -= fac*Ni_p[i]*fz[p];
787: }
788: }
789: }
791: static void FormContinuityRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar hc[])
792: {
793: PetscInt ngp;
794: PetscScalar gp_xi[GAUSS_POINTS][NSD];
795: PetscScalar gp_weight[GAUSS_POINTS];
796: PetscInt p,i;
797: PetscScalar Ni_p[NODES_PER_EL];
798: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
799: PetscScalar J_p,fac;
801: /* define quadrature rule */
802: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
804: /* evaluate integral */
805: for (p = 0; p < ngp; p++) {
806: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
807: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
808: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
809: fac = gp_weight[p]*J_p;
811: for (i = 0; i < NODES_PER_EL; i++) {
812: Fe[i] -= fac*Ni_p[i]*hc[p];
813: }
814: }
815: }
817: #define _ZERO_ROWCOL_i(A,i) { \
818: PetscInt KK; \
819: PetscScalar tmp = A[24*(i)+(i)]; \
820: for(KK=0;KK<24;KK++){A[24*(i)+KK]=0.0;} \
821: for(KK=0;KK<24;KK++){A[24*KK+(i)]=0.0;} \
822: A[24*(i)+(i)] = tmp;} \
824: #define _ZERO_ROW_i(A,i) { \
825: PetscInt KK; \
826: for(KK=0;KK<8;KK++){A[8*(i)+KK]=0.0;}}
828: #define _ZERO_COL_i(A,i) { \
829: PetscInt KK; \
830: for(KK=0;KK<8;KK++){A[24*KK+(i)]=0.0;}}
834: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,CellProperties cell_properties)
835: {
836: DM cda;
837: Vec coords;
838: DMDACoor3d ***_coords;
839: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
840: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
841: PetscInt sex,sey,sez,mx,my,mz;
842: PetscInt ei,ej,ek;
843: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
844: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
845: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
846: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
847: PetscScalar el_coords[NODES_PER_EL*NSD];
848: GaussPointCoefficients *props;
849: PetscScalar *prop_eta;
850: PetscInt n,M,N,P;
851: PetscLogDouble t0,t1;
852: PetscErrorCode ierr;
856: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
857: /* setup for coords */
858: DMDAGetCoordinateDA(stokes_da,&cda);
859: DMDAGetGhostedCoordinates(stokes_da,&coords);
860: DMDAVecGetArray(cda,coords,&_coords);
862: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
863: PetscGetTime(&t0);
864: for (ek = sez; ek < sez+mz; ek++) {
865: for (ej = sey; ej < sey+my; ej++) {
866: for (ei = sex; ei < sex+mx; ei++) {
867: /* get coords for the element */
868: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
869: /* get cell properties */
870: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
871: /* get coefficients for the element */
872: prop_eta = props->eta;
874: /* initialise element stiffness matrix */
875: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
876: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
877: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
878: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
880: /* form element stiffness matrix */
881: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
882: FormGradientOperatorQ13D(Ge,el_coords);
883: /*#ifdef ASSEMBLE_LOWER_TRIANGULAR*/
884: FormDivergenceOperatorQ13D(De,el_coords);
885: /*#endif*/
886: FormStabilisationOperatorQ13D(Ce,el_coords,prop_eta);
888: /* insert element matrix into global matrix */
889: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
891: for (n=0; n<NODES_PER_EL; n++) {
892: if ( (u_eqn[3*n ].i == 0) || (u_eqn[3*n ].i == M-1) ) {
893: _ZERO_ROWCOL_i(Ae,3*n);
894: _ZERO_ROW_i (Ge,3*n);
895: _ZERO_COL_i (De,3*n);
896: }
898: if ( (u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1) ) {
899: _ZERO_ROWCOL_i(Ae,3*n+1);
900: _ZERO_ROW_i (Ge,3*n+1);
901: _ZERO_COL_i (De,3*n+1);
902: }
904: if ( (u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1) ) {
905: _ZERO_ROWCOL_i(Ae,3*n+2);
906: _ZERO_ROW_i (Ge,3*n+2);
907: _ZERO_COL_i (De,3*n+2);
908: }
909: }
910: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
911: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
912: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
913: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
914: }
915: }
916: }
917: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
918: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
920: DMDAVecRestoreArray(cda,coords,&_coords);
922: PetscGetTime(&t1);
923: return(0);
924: }
928: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,CellProperties cell_properties)
929: {
930: DM cda;
931: Vec coords;
932: DMDACoor3d ***_coords;
933: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
934: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
935: PetscInt sex,sey,sez,mx,my,mz;
936: PetscInt ei,ej,ek;
937: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
938: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
939: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
940: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
941: PetscScalar el_coords[NODES_PER_EL*NSD];
942: GaussPointCoefficients *props;
943: PetscScalar *prop_eta;
944: PetscInt n,M,N,P;
945: PetscErrorCode ierr;
949: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
950: /* setup for coords */
951: DMDAGetCoordinateDA(stokes_da,&cda);
952: DMDAGetGhostedCoordinates(stokes_da,&coords);
953: DMDAVecGetArray(cda,coords,&_coords);
955: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
956: for (ek = sez; ek < sez+mz; ek++) {
957: for (ej = sey; ej < sey+my; ej++) {
958: for (ei = sex; ei < sex+mx; ei++) {
959: /* get coords for the element */
960: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
961: /* get cell properties */
962: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
963: /* get coefficients for the element */
964: prop_eta = props->eta;
966: /* initialise element stiffness matrix */
967: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
968: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
969: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
970: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
972: /* form element stiffness matrix */
973: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
974: FormGradientOperatorQ13D(Ge,el_coords);
975: /* FormDivergenceOperatorQ13D(De,el_coords); */
976: FormScaledMassMatrixOperatorQ13D(Ce,el_coords,prop_eta);
978: /* insert element matrix into global matrix */
979: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
981: for (n=0; n<NODES_PER_EL; n++) {
982: if ( (u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1) ) {
983: _ZERO_ROWCOL_i(Ae,3*n);
984: _ZERO_ROW_i (Ge,3*n);
985: _ZERO_COL_i (De,3*n);
986: }
988: if ( (u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1) ) {
989: _ZERO_ROWCOL_i(Ae,3*n+1);
990: _ZERO_ROW_i (Ge,3*n+1);
991: _ZERO_COL_i (De,3*n+1);
992: }
994: if ( (u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1) ) {
995: _ZERO_ROWCOL_i(Ae,3*n+2);
996: _ZERO_ROW_i (Ge,3*n+2);
997: _ZERO_COL_i (De,3*n+2);
998: }
999: }
1000: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
1001: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
1002: /*MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);*/
1003: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
1004: }
1005: }
1006: }
1007: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
1008: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
1010: DMDAVecRestoreArray(cda,coords,&_coords);
1011: return(0);
1012: }
1016: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,CellProperties cell_properties)
1017: {
1018: DM cda;
1019: Vec coords;
1020: DMDACoor3d ***_coords;
1021: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
1022: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
1023: PetscInt sex,sey,sez,mx,my,mz;
1024: PetscInt ei,ej,ek;
1025: PetscScalar Fe[NODES_PER_EL*U_DOFS];
1026: PetscScalar He[NODES_PER_EL*P_DOFS];
1027: PetscScalar el_coords[NODES_PER_EL*NSD];
1028: GaussPointCoefficients *props;
1029: PetscScalar *prop_fx,*prop_fy,*prop_fz,*prop_hc;
1030: Vec local_F;
1031: StokesDOF ***ff;
1032: PetscInt n,M,N,P;
1033: PetscErrorCode ierr;
1037: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
1038: /* setup for coords */
1039: DMDAGetCoordinateDA(stokes_da,&cda);
1040: DMDAGetGhostedCoordinates(stokes_da,&coords);
1041: DMDAVecGetArray(cda,coords,&_coords);
1043: /* get acces to the vector */
1044: DMGetLocalVector(stokes_da,&local_F);
1045: VecZeroEntries(local_F);
1046: DMDAVecGetArray(stokes_da,local_F,&ff);
1048: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1049: for (ek = sez; ek < sez+mz; ek++) {
1050: for (ej = sey; ej < sey+my; ej++) {
1051: for (ei = sex; ei < sex+mx; ei++) {
1052: /* get coords for the element */
1053: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1054: /* get cell properties */
1055: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
1056: /* get coefficients for the element */
1057: prop_fx = props->fx;
1058: prop_fy = props->fy;
1059: prop_fz = props->fz;
1060: prop_hc = props->hc;
1062: /* initialise element stiffness matrix */
1063: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
1064: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
1066: /* form element stiffness matrix */
1067: FormMomentumRhsQ13D(Fe,el_coords,prop_fx,prop_fy,prop_fz);
1068: FormContinuityRhsQ13D(He,el_coords,prop_hc);
1070: /* insert element matrix into global matrix */
1071: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
1073: for (n=0; n<NODES_PER_EL; n++) {
1074: if ( (u_eqn[3*n ].i == 0) || (u_eqn[3*n ].i == M-1) ) {
1075: Fe[3*n ] = 0.0;
1076: }
1078: if ( (u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1) ) {
1079: Fe[3*n+1] = 0.0;
1080: }
1082: if ( (u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1) ) {
1083: Fe[3*n+2] = 0.0;
1084: }
1085: }
1087: DMDASetValuesLocalStencil3D_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
1088: }
1089: }
1090: }
1091: DMDAVecRestoreArray(stokes_da,local_F,&ff);
1092: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
1093: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
1094: DMRestoreLocalVector(stokes_da,&local_F);
1096: DMDAVecRestoreArray(cda,coords,&_coords);
1097: return(0);
1098: }
1100: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta,PetscReal *MX,PetscReal *MY,PetscReal *MZ)
1101: {
1102: *theta = 0.0;
1103: *MX = 2.0 * M_PI;
1104: *MY = 2.0 * M_PI;
1105: *MZ = 2.0 * M_PI;
1106: }
1107: static void evaluate_MS_FrankKamentski(PetscReal pos[],PetscReal v[],PetscReal *p,PetscReal *eta,PetscReal Fm[],PetscReal *Fc)
1108: {
1109: PetscReal x,y,z;
1110: PetscReal theta,MX,MY,MZ;
1112: evaluate_MS_FrankKamentski_constants(&theta,&MX,&MY,&MZ);
1113: x = pos[0];
1114: y = pos[1];
1115: z = pos[2];
1116: if (v) {
1117: /*
1118: v[0] = pow(z,3)*exp(y)*sin(M_PI*x);
1119: v[1] = z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y);
1120: v[2] = -(pow(x,3) + pow(y,3))*sin(2.0*M_PI*z);
1121: */
1122: /*
1123: v[0] = sin(M_PI*x);
1124: v[1] = sin(M_PI*y);
1125: v[2] = sin(M_PI*z);
1126: */
1127: v[0] = pow(z,3)*exp(y)*sin(M_PI*x);
1128: v[1] = z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y);
1129: v[2] = pow(z,2)*(cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 - M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2) - M_PI*pow(z,4)*cos(M_PI*x)*exp(y)/4 ;
1130: }
1131: if (p) {
1132: *p = pow(x,2) + pow(y,2) + pow(z,2);
1133: }
1134: if (eta) {
1135: /**eta = exp(-theta*(1.0 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)));*/
1136: *eta = 1.0;
1137: }
1138: if (Fm) {
1139: /*
1140: Fm[0] = -2*x - 2.0*pow(M_PI,2)*pow(z,3)*exp(y)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(M_PI*x) - 0.2*MZ*theta*(-1.5*pow(x,2)*sin(2.0*M_PI*z) + 1.5*pow(z,2)*exp(y)*sin(M_PI*x))*cos(MX*x)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MY*y)*sin(MZ*z) - 0.2*M_PI*MX*theta*pow(z,3)*cos(M_PI*x)*cos(MZ*z)*exp(y)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MX*x)*sin(MY*y) + 2.0*(3.0*z*exp(y)*sin(M_PI*x) - 3.0*M_PI*pow(x,2)*cos(2.0*M_PI*z))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*(0.5*pow(z,3)*exp(y)*sin(M_PI*x) + M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x) - 1.0*z*pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*theta*(1 + 0.1*MY*cos(MX*x)*cos(MY*y)*cos(MZ*z))*(0.5*pow(z,3)*exp(y)*sin(M_PI*x) - 1.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) ;
1141: Fm[1] = -2*y - 0.2*MX*theta*(0.5*pow(z,3)*exp(y)*sin(M_PI*x) - 1.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x))*cos(MZ*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MX*x)*sin(MY*y) - 0.2*MZ*theta*(-1.5*pow(y,2)*sin(2.0*M_PI*z) + 0.5*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y))*cos(MX*x)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MY*y)*sin(MZ*z) + 2.0*(-2.0*z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 0.5*M_PI*pow(z,3)*cos(M_PI*x)*exp(y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*(z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 2*M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*theta*(1 + 0.1*MY*cos(MX*x)*cos(MY*y)*cos(MZ*z))*(-z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) - 6.0*M_PI*pow(y,2)*cos(2.0*M_PI*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)));
1142: Fm[2] = -2*z + 8.0*pow(M_PI,2)*(pow(x,3) + pow(y,3))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(2.0*M_PI*z) - 0.2*MX*theta*(-1.5*pow(x,2)*sin(2.0*M_PI*z) + 1.5*pow(z,2)*exp(y)*sin(M_PI*x))*cos(MZ*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MX*x)*sin(MY*y) + 0.4*M_PI*MZ*theta*(pow(x,3) + pow(y,3))*cos(MX*x)*cos(2.0*M_PI*z)*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)))*sin(MY*y)*sin(MZ*z) + 2.0*(-3.0*x*sin(2.0*M_PI*z) + 1.5*M_PI*pow(z,2)*cos(M_PI*x)*exp(y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*(-3.0*y*sin(2.0*M_PI*z) - 0.5*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 0.5*M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) + 2.0*theta*(1 + 0.1*MY*cos(MX*x)*cos(MY*y)*cos(MZ*z))*(-1.5*pow(y,2)*sin(2.0*M_PI*z) + 0.5*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y))*exp(-theta*(1 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y))) ;
1143: */
1144: /*
1145: Fm[0]=-2*x - 2.0*pow(M_PI,2)*sin(M_PI*x);
1146: Fm[1]=-2*y - 2.0*pow(M_PI,2)*sin(M_PI*y);
1147: Fm[2]=-2*z - 2.0*pow(M_PI,2)*sin(M_PI*z);
1148: */
1149: /*
1150: Fm[0] = -2*x + pow(z,3)*exp(y)*sin(M_PI*x) + 6.0*z*exp(y)*sin(M_PI*x) - 6.0*M_PI*pow(x,2)*cos(2.0*M_PI*z) - 2.0*pow(M_PI,2)*pow(z,3)*exp(y)*sin(M_PI*x) + 2.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x) - 2.0*z*pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x) ;
1151: Fm[1] = -2*y - 6.0*M_PI*pow(y,2)*cos(2.0*M_PI*z) + 2.0*z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 6.0*z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*pow(z,3)*cos(M_PI*x)*exp(y) - 4.0*M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y);
1152: Fm[2] = -2*z - 6.0*x*sin(2.0*M_PI*z) - 6.0*y*sin(2.0*M_PI*z) - cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 8.0*pow(M_PI,2)*(pow(x,3) + pow(y,3))*sin(2.0*M_PI*z) + 3.0*M_PI*pow(z,2)*cos(M_PI*x)*exp(y) + M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;
1153: */
1154: Fm[0] = -2*x + 2*z*(pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x) - 1.0*M_PI*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x)) + pow(z,3)*exp(y)*sin(M_PI*x) + 6.0*z*exp(y)*sin(M_PI*x) - 1.0*pow(M_PI,2)*pow(z,3)*exp(y)*sin(M_PI*x) + 2.0*M_PI*z*exp(-y)*sin(M_PI*y)*sin(2.0*M_PI*x) - 2.0*z*pow(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x);
1155: Fm[1] = -2*y + 2*z*(-cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)) + 2.0*z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 6.0*z*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) - 4.0*M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;
1156: Fm[2] = -2*z + pow(z,2)*(-2.0*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 2.0*pow(M_PI,3)*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)) + pow(z,2)*(cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 - 3*pow(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + pow(M_PI,3)*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2 - 3*M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)/2) + 1.0*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 0.25*pow(M_PI,3)*pow(z,4)*cos(M_PI*x)*exp(y) - 0.25*M_PI*pow(z,4)*cos(M_PI*x)*exp(y) - 3.0*M_PI*pow(z,2)*cos(M_PI*x)*exp(y) - 1.0*M_PI*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;
1157: }
1158: if (Fc) {
1159: /**Fc = -2.0*M_PI*(pow(x,3) + pow(y,3))*cos(2.0*M_PI*z) - z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*pow(z,3)*cos(M_PI*x)*exp(y) + M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y) ;*/
1160: /**Fc = M_PI*cos(M_PI*x) + M_PI*cos(M_PI*y) + M_PI*cos(M_PI*z);*/
1161: /**Fc = -2.0*M_PI*(pow(x,3) + pow(y,3))*cos(2.0*M_PI*z) - z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + M_PI*pow(z,3)*cos(M_PI*x)*exp(y) + M_PI*z*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y);*/
1162: *Fc = 0.0;
1163: }
1164: }
1168: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx,PetscInt my,PetscInt mz,DM *_da,Vec *_X)
1169: {
1170: DM da,cda;
1171: Vec X,local_X;
1172: StokesDOF ***_stokes;
1173: Vec coords;
1174: DMDACoor3d ***_coords;
1175: PetscInt si,sj,sk,ei,ej,ek,i,j,k;
1179: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1180: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,4,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&da);
1181: DMDASetFieldName(da,0,"anlytic_Vx");
1182: DMDASetFieldName(da,1,"anlytic_Vy");
1183: DMDASetFieldName(da,2,"anlytic_Vz");
1184: DMDASetFieldName(da,3,"analytic_P");
1186: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
1188: DMDAGetGhostedCoordinates(da,&coords);
1189: DMDAGetCoordinateDA(da,&cda);
1190: DMDAVecGetArray(cda,coords,&_coords);
1192: DMCreateGlobalVector(da,&X);
1193: DMCreateLocalVector(da,&local_X);
1194: DMDAVecGetArray(da,local_X,&_stokes);
1196: DMDAGetGhostCorners(da,&si,&sj,&sk,&ei,&ej,&ek);
1197: for (k = sk; k < sk+ek; k++) {
1198: for (j = sj; j < sj+ej; j++) {
1199: for (i = si; i < si+ei; i++) {
1200: PetscReal pos[NSD],pressure,vel[NSD];
1202: pos[0] = PetscRealPart(_coords[k][j][i].x);
1203: pos[1] = PetscRealPart(_coords[k][j][i].y);
1204: pos[2] = PetscRealPart(_coords[k][j][i].z);
1206: evaluate_MS_FrankKamentski(pos,vel,&pressure,PETSC_NULL,PETSC_NULL,PETSC_NULL);
1208: _stokes[k][j][i].u_dof = vel[0];
1209: _stokes[k][j][i].v_dof = vel[1];
1210: _stokes[k][j][i].w_dof = vel[2];
1211: _stokes[k][j][i].p_dof = pressure;
1212: }
1213: }
1214: }
1215: DMDAVecRestoreArray(da,local_X,&_stokes);
1216: DMDAVecRestoreArray(cda,coords,&_coords);
1218: DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1219: DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);
1221: VecDestroy(&local_X);
1223: *_da = da;
1224: *_X = X;
1225: return(0);
1226: }
1230: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da,Vec X,Vec X_analytic)
1231: {
1232: DM cda;
1233: Vec coords,X_analytic_local,X_local;
1234: DMDACoor3d ***_coords;
1235: PetscInt sex,sey,sez,mx,my,mz;
1236: PetscInt ei,ej,ek;
1237: PetscScalar el_coords[NODES_PER_EL*NSD];
1238: StokesDOF ***stokes_analytic,***stokes;
1239: StokesDOF stokes_analytic_e[NODES_PER_EL],stokes_e[NODES_PER_EL];
1241: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1242: PetscScalar Ni_p[NODES_PER_EL];
1243: PetscInt ngp;
1244: PetscScalar gp_xi[GAUSS_POINTS][NSD];
1245: PetscScalar gp_weight[GAUSS_POINTS];
1246: PetscInt p,i;
1247: PetscScalar J_p,fac;
1248: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1249: PetscScalar tint_p_ms,tint_p,int_p_ms,int_p;
1250: PetscInt M;
1251: PetscReal xymin[NSD],xymax[NSD];
1255: /* define quadrature rule */
1256: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1258: /* setup for coords */
1259: DMDAGetCoordinateDA(stokes_da,&cda);
1260: DMDAGetGhostedCoordinates(stokes_da,&coords);
1261: DMDAVecGetArray(cda,coords,&_coords);
1263: /* setup for analytic */
1264: DMCreateLocalVector(stokes_da,&X_analytic_local);
1265: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1266: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1267: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1269: /* setup for solution */
1270: DMCreateLocalVector(stokes_da,&X_local);
1271: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1272: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1273: DMDAVecGetArray(stokes_da,X_local,&stokes);
1275: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1276: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1278: h = (xymax[0]-xymin[0])/((PetscReal)(M-1));
1280: tp_L2 = tu_L2 = tu_H1 = 0.0;
1281: tint_p_ms = tint_p = 0.0;
1283: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1285: for (ek = sez; ek < sez+mz; ek++) {
1286: for (ej = sey; ej < sey+my; ej++) {
1287: for (ei = sex; ei < sex+mx; ei++) {
1288: /* get coords for the element */
1289: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1290: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1291: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1293: /* evaluate integral */
1294: for (p = 0; p < ngp; p++) {
1295: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1296: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1297: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1298: fac = gp_weight[p]*J_p;
1300: for (i = 0; i < NODES_PER_EL; i++) {
1301: tint_p_ms = tint_p_ms+fac*Ni_p[i]*stokes_analytic_e[i].p_dof;
1302: tint_p = tint_p +fac*Ni_p[i]*stokes_e[i].p_dof;
1303: }
1304: }
1305: }
1306: }
1307: }
1308: MPI_Allreduce(&tint_p_ms,&int_p_ms,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1309: MPI_Allreduce(&tint_p,&int_p,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1311: PetscPrintf(PETSC_COMM_WORLD,"\\int P dv %1.4e (h) %1.4e (ms)\n",PetscRealPart(int_p),PetscRealPart(int_p_ms));
1313: /* remove mine and add manufacture one */
1314: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1315: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1317: {
1318: PetscInt k,L,dof;
1319: PetscScalar *fields;
1320: DMDAGetInfo(stokes_da,0, 0,0,0, 0,0,0, &dof,0,0,0,0,0);
1322: VecGetLocalSize(X_local,&L);
1323: VecGetArray(X_local,&fields);
1324: for (k=0; k<L/dof; k++) {
1325: fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1326: }
1327: VecRestoreArray(X_local,&fields);
1329: VecGetLocalSize(X,&L);
1330: VecGetArray(X,&fields);
1331: for (k=0; k<L/dof; k++) {
1332: fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1333: }
1334: VecRestoreArray(X,&fields);
1335: }
1337: DMDAVecGetArray(stokes_da,X_local,&stokes);
1338: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1340: for (ek = sez; ek < sez+mz; ek++) {
1341: for (ej = sey; ej < sey+my; ej++) {
1342: for (ei = sex; ei < sex+mx; ei++) {
1343: /* get coords for the element */
1344: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1345: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1346: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1348: /* evaluate integral */
1349: p_e_L2 = 0.0;
1350: u_e_L2 = 0.0;
1351: u_e_H1 = 0.0;
1352: for (p = 0; p < ngp; p++) {
1353: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1354: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1355: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1356: fac = gp_weight[p]*J_p;
1358: for (i = 0; i < NODES_PER_EL; i++) {
1359: PetscScalar u_error,v_error,w_error;
1361: 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);
1363: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1364: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1365: w_error = stokes_e[i].w_dof-stokes_analytic_e[i].w_dof;
1366: /*
1367: if (p==0) {
1368: printf("p=0: %d %d %d %1.4e,%1.4e,%1.4e \n", ei,ej,ek,u_error,v_error,w_error);
1369: }
1370: */
1371: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error+w_error*w_error);
1373: u_e_H1 = u_e_H1+fac*( GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1374: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error
1375: +GNx_p[2][i]*u_error*GNx_p[2][i]*u_error
1376: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1377: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error
1378: +GNx_p[2][i]*v_error*GNx_p[2][i]*v_error
1379: +GNx_p[0][i]*w_error*GNx_p[0][i]*w_error /* dw/dx */
1380: +GNx_p[1][i]*w_error*GNx_p[1][i]*w_error
1381: +GNx_p[2][i]*w_error*GNx_p[2][i]*w_error );
1382: }
1383: }
1385: tp_L2 += p_e_L2;
1386: tu_L2 += u_e_L2;
1387: tu_H1 += u_e_H1;
1388: }
1389: }
1390: }
1391: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1392: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1393: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1394: p_L2 = PetscSqrtScalar(p_L2);
1395: u_L2 = PetscSqrtScalar(u_L2);
1396: u_H1 = PetscSqrtScalar(u_H1);
1398: PetscPrintf(PETSC_COMM_WORLD,"%1.4e %1.4e %1.4e %1.4e \n",PetscRealPart(h),PetscRealPart(p_L2),PetscRealPart(u_L2),PetscRealPart(u_H1));
1401: DMDAVecRestoreArray(cda,coords,&_coords);
1403: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1404: VecDestroy(&X_analytic_local);
1405: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1406: VecDestroy(&X_local);
1407: return(0);
1408: }
1412: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da,Vec FIELD,const char file_prefix[])
1413: {
1414: char vtk_filename[PETSC_MAX_PATH_LEN];
1415: PetscMPIInt rank;
1416: MPI_Comm comm;
1417: FILE *vtk_fp = PETSC_NULL;
1418: PetscInt si,sj,sk,nx,ny,nz,i;
1419: PetscInt f,n_fields,N;
1420: DM cda;
1421: Vec coords;
1422: Vec l_FIELD;
1423: PetscScalar *_L_FIELD;
1424: PetscInt memory_offset;
1425: PetscScalar *buffer;
1426: PetscLogDouble t0,t1;
1430: PetscGetTime(&t0);
1432: /* create file name */
1433: PetscObjectGetComm((PetscObject)da,&comm);
1434: MPI_Comm_rank(comm,&rank);
1435: PetscSNPrintf(vtk_filename,sizeof vtk_filename,"subdomain-%s-p%1.4d.vts",file_prefix,rank);
1437: /* open file and write header */
1438: vtk_fp = fopen(vtk_filename,"w");
1439: if (!vtk_fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename);
1441: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");
1443: /* coords */
1444: DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1445: N = nx * ny * nz;
1447: #ifdef PETSC_WORDS_BIGENDIAN
1448: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1449: #else
1450: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1451: #endif
1452: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <StructuredGrid WholeExtent=\"%D %D %D %D %D %D\">\n",si,si+nx-1,sj,sj+ny-1,sk,sk+nz-1);
1453: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <Piece Extent=\"%D %D %D %D %D %D\">\n",si,si+nx-1,sj,sj+ny-1,sk,sk+nz-1);
1455: memory_offset = 0;
1457: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <CellData></CellData>\n");
1459: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <Points>\n");
1461: /* copy coordinates */
1462: DMDAGetCoordinateDA(da,&cda);
1463: DMDAGetGhostedCoordinates(da,&coords);
1464: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%d\" />\n",memory_offset);
1465: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N*3;
1467: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </Points>\n");
1469: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PointData Scalars=\" ");
1470: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_fields,0,0,0,0,0);
1471: for (f=0; f<n_fields; f++) {
1472: const char *field_name;
1473: DMDAGetFieldName(da,f,&field_name);
1474: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"%s ",field_name);
1475: }
1476: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\">\n");
1478: for (f=0; f<n_fields; f++) {
1479: const char *field_name;
1481: DMDAGetFieldName(da,f,&field_name);
1482: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <DataArray type=\"Float64\" Name=\"%s\" format=\"appended\" offset=\"%d\"/>\n", field_name,memory_offset);
1483: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N;
1484: }
1486: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PointData>\n");
1487: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </Piece>\n");
1488: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </StructuredGrid>\n");
1490: PetscMalloc(sizeof(PetscScalar)*N,&buffer);
1491: DMGetLocalVector(da,&l_FIELD);
1492: DMGlobalToLocalBegin(da, FIELD,INSERT_VALUES,l_FIELD);
1493: DMGlobalToLocalEnd(da,FIELD,INSERT_VALUES,l_FIELD);
1494: VecGetArray(l_FIELD,&_L_FIELD);
1496: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <AppendedData encoding=\"raw\">\n");
1497: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"_");
1499: /* write coordinates */
1500: {
1501: int length = sizeof(PetscScalar)*N*3;
1502: PetscScalar *allcoords;
1504: fwrite(&length,sizeof(int),1,vtk_fp);
1505: VecGetArray(coords,&allcoords);
1506: fwrite(allcoords,sizeof(PetscScalar),3*N,vtk_fp);
1507: VecRestoreArray(coords,&allcoords);
1508: }
1509: /* write fields */
1510: for (f=0; f<n_fields; f++) {
1511: int length = sizeof(PetscScalar)*N;
1512: fwrite(&length,sizeof(int),1,vtk_fp);
1513: /* load */
1514: for (i=0; i<N; i++) {
1515: buffer[i] = _L_FIELD[n_fields*i + f];
1516: }
1518: /* write */
1519: fwrite(buffer,sizeof(PetscScalar),N,vtk_fp);
1520: }
1521: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\n </AppendedData>\n");
1523: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");
1525: PetscFree(buffer);
1526: VecRestoreArray(l_FIELD,&_L_FIELD);
1527: DMRestoreLocalVector(da,&l_FIELD);
1529: if (vtk_fp) {
1530: fclose(vtk_fp);
1531: vtk_fp = NULL;
1532: }
1534: PetscGetTime(&t1);
1535: return(0);
1536: }
1540: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp,PetscInt indent_level,DM da,const char local_file_prefix[])
1541: {
1542: PetscMPIInt nproc,rank;
1543: MPI_Comm comm;
1544: const PetscInt *lx,*ly,*lz;
1545: PetscInt M,N,P,pM,pN,pP,sum,*olx,*oly,*olz;
1546: PetscInt *osx,*osy,*osz,*oex,*oey,*oez;
1547: PetscInt i,j,k,II,stencil;
1551: /* create file name */
1552: PetscObjectGetComm((PetscObject)da,&comm);
1553: MPI_Comm_size(comm,&nproc);
1554: MPI_Comm_rank(comm,&rank);
1556: DMDAGetInfo(da,0,&M,&N,&P,&pM,&pN,&pP,0,&stencil,0,0,0,0);
1557: DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
1559: /* generate start,end list */
1560: PetscMalloc(sizeof(PetscInt)*(pM+1),&olx);
1561: PetscMalloc(sizeof(PetscInt)*(pN+1),&oly);
1562: PetscMalloc(sizeof(PetscInt)*(pP+1),&olz);
1563: sum = 0;
1564: for (i=0; i<pM; i++) {
1565: olx[i] = sum;
1566: sum = sum + lx[i];
1567: } olx[pM] = sum;
1568: sum = 0;
1569: for (i=0; i<pN; i++) {
1570: oly[i] = sum;
1571: sum = sum + ly[i];
1572: } oly[pN] = sum;
1573: sum = 0;
1574: for (i=0; i<pP; i++) {
1575: olz[i] = sum;
1576: sum = sum + lz[i];
1577: } olz[pP] = sum;
1578: PetscMalloc(sizeof(PetscInt)*(pM),&osx);
1579: PetscMalloc(sizeof(PetscInt)*(pN),&osy);
1580: PetscMalloc(sizeof(PetscInt)*(pP),&osz);
1581: PetscMalloc(sizeof(PetscInt)*(pM),&oex);
1582: PetscMalloc(sizeof(PetscInt)*(pN),&oey);
1583: PetscMalloc(sizeof(PetscInt)*(pP),&oez);
1584: for (i=0; i<pM; i++) {
1585: osx[i] = olx[i] - stencil;
1586: oex[i] = olx[i] + lx[i] + stencil;
1587: if (osx[i]<0)osx[i]=0;
1588: if (oex[i]>M)oex[i]=M;
1589: }
1591: for (i=0; i<pN; i++) {
1592: osy[i] = oly[i] - stencil;
1593: oey[i] = oly[i] + ly[i] + stencil;
1594: if (osy[i]<0)osy[i]=0;
1595: if (oey[i]>M)oey[i]=N;
1596: }
1597: for (i=0; i<pP; i++) {
1598: osz[i] = olz[i] - stencil;
1599: oez[i] = olz[i] + lz[i] + stencil;
1600: if (osz[i]<0)osz[i]=0;
1601: if (oez[i]>P)oez[i]=P;
1602: }
1604: for (k=0; k<pP; k++) {
1605: for (j=0; j<pN; j++) {
1606: for (i=0; i<pM; i++) {
1607: char name[PETSC_MAX_PATH_LEN];
1608: PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
1609: PetscSNPrintf(name,sizeof name,"subdomain-%s-p%1.4d.vts",local_file_prefix,procid);
1610: for (II=0; II<indent_level; II++) {
1611: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," ");
1612: }
1613: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<Piece Extent=\"%d %d %d %d %d %d\" Source=\"%s\"/>\n",
1614: osx[i],oex[i]-1,
1615: osy[j],oey[j]-1,
1616: osz[k],oez[k]-1,name);
1617: }
1618: }
1619: }
1620: PetscFree(olx);
1621: PetscFree(oly);
1622: PetscFree(olz);
1623: PetscFree(osx);
1624: PetscFree(osy);
1625: PetscFree(osz);
1626: PetscFree(oex);
1627: PetscFree(oey);
1628: PetscFree(oez);
1629: return(0);
1630: }
1634: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da,const char file_prefix[],const char local_file_prefix[])
1635: {
1636: MPI_Comm comm;
1637: PetscMPIInt nproc,rank;
1638: char vtk_filename[PETSC_MAX_PATH_LEN];
1639: FILE *vtk_fp = PETSC_NULL;
1640: PetscInt M,N,P,si,sj,sk,nx,ny,nz;
1641: PetscInt i,dofs;
1645: /* only master generates this file */
1646: PetscObjectGetComm((PetscObject)da,&comm);
1647: MPI_Comm_size(comm,&nproc);
1648: MPI_Comm_rank(comm,&rank);
1650: if (rank != 0) { return(0); }
1652: /* create file name */
1653: PetscSNPrintf(vtk_filename,sizeof vtk_filename,"%s.pvts",file_prefix);
1654: vtk_fp = fopen(vtk_filename,"w");
1655: if (!vtk_fp) { SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename); }
1657: /* (VTK) generate pvts header */
1658: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");
1660: #ifdef PETSC_WORDS_BIGENDIAN
1661: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1662: #else
1663: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1664: #endif
1666: /* define size of the nodal mesh based on the cell DM */
1667: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,&dofs,0,0,0,0,0);
1668: DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1669: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PStructuredGrid GhostLevel=\"1\" WholeExtent=\"%d %d %d %d %d %d\">\n",0,M-1,0,N-1,0,P-1); /* note overlap = 1 for Q1 */
1671: /* DUMP THE CELL REFERENCES */
1672: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PCellData>\n");
1673: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PCellData>\n");
1675: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PPoints>\n");
1676: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"%d\"/>\n",NSD);
1677: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PPoints>\n");
1679: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PPointData>\n");
1680: for (i=0; i<dofs; i++) {
1681: const char *fieldname;
1682: DMDAGetFieldName(da,i,&fieldname);
1683: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n",fieldname);
1684: }
1685: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PPointData>\n");
1687: /* write out the parallel information */
1688: DAViewVTK_write_PieceExtend(vtk_fp,2,da,local_file_prefix);
1690: /* close the file */
1691: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PStructuredGrid>\n");
1692: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");
1694: if (vtk_fp){
1695: fclose(vtk_fp);
1696: vtk_fp = NULL;
1697: }
1698: return(0);
1699: }
1703: PetscErrorCode DAView3DPVTS(DM da, Vec x,const char NAME[])
1704: {
1705: char vts_filename[PETSC_MAX_PATH_LEN];
1706: char pvts_filename[PETSC_MAX_PATH_LEN];
1710: PetscSNPrintf(vts_filename,sizeof vts_filename,"%s-mesh",NAME);
1711: DAView_3DVTK_StructuredGrid_appended(da,x,vts_filename);
1713: PetscSNPrintf(pvts_filename,sizeof pvts_filename,"%s-mesh",NAME);
1714: DAView_3DVTK_PStructuredGrid(da,pvts_filename,vts_filename);
1715: return(0);
1716: }
1720: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
1721: {
1723: PetscReal norms[4];
1724: Vec Br,v,w;
1725: Mat A;
1728: KSPGetOperators(ksp,&A,0,0);
1729: MatGetVecs(A,&w,&v);
1731: KSPBuildResidual(ksp,v,w,&Br);
1733: VecStrideNorm(Br,0,NORM_2,&norms[0]);
1734: VecStrideNorm(Br,1,NORM_2,&norms[1]);
1735: VecStrideNorm(Br,2,NORM_2,&norms[2]);
1736: VecStrideNorm(Br,3,NORM_2,&norms[3]);
1738: VecDestroy(&v);
1739: VecDestroy(&w);
1741: PetscPrintf(PETSC_COMM_WORLD,"%3D KSP Component U,V,W,P residual norm [ %1.12e, %1.12e, %1.12e, %1.12e ]\n",n,norms[0],norms[1],norms[2],norms[3]);
1742: return(0);
1743: }
1747: static PetscErrorCode PCMGSetupViaCoarsen(PC pc,DM da_fine)
1748: {
1749: PetscInt nlevels,k,PETSC_UNUSED finest;
1750: DM *da_list,*daclist;
1751: Mat R;
1755: nlevels = 1;
1756: PetscOptionsGetInt(PETSC_NULL,"-levels",&nlevels,0);
1758: PetscMalloc(sizeof(DM)*nlevels,&da_list);
1759: for (k=0; k<nlevels; k++) { da_list[k] = PETSC_NULL; }
1760: PetscMalloc(sizeof(DM)*nlevels,&daclist);
1761: for (k=0; k<nlevels; k++) { daclist[k] = PETSC_NULL; }
1763: /* finest grid is nlevels - 1 */
1764: finest = nlevels - 1;
1765: daclist[0] = da_fine;
1766: PetscObjectReference((PetscObject)da_fine);
1767: DMCoarsenHierarchy(da_fine,nlevels-1,&daclist[1]);
1768: for (k=0; k<nlevels; k++) {
1769: da_list[k] = daclist[nlevels-1-k];
1770: DMDASetUniformCoordinates(da_list[k],0.0,1.0,0.0,1.0,0.0,1.0);
1771: }
1773: PCMGSetLevels(pc,nlevels,PETSC_NULL);
1774: PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
1775: PCMGSetGalerkin(pc,PETSC_TRUE);
1777: for (k=1; k<nlevels; k++){
1778: DMCreateInterpolation(da_list[k-1],da_list[k],&R,PETSC_NULL);
1779: PCMGSetInterpolation(pc,k,R);
1780: MatDestroy(&R);
1781: }
1783: /* tidy up */
1784: for (k=0; k<nlevels; k++) {
1785: DMDestroy(&da_list[k]);
1786: }
1787: PetscFree(da_list);
1788: PetscFree(daclist);
1789: return(0);
1790: }
1794: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx,PetscInt my,PetscInt mz)
1795: {
1796: DM da_Stokes;
1797: PetscInt u_dof,p_dof,dof,stencil_width;
1798: Mat A,B;
1799: PetscInt mxl,myl,mzl;
1800: DM vel_cda;
1801: Vec vel_coords;
1802: PetscInt p;
1803: Vec f,X;
1804: DMDACoor3d ***_vel_coords;
1805: PetscInt its;
1806: KSP ksp_S;
1807: PetscInt model_definition = 0;
1808: PetscInt ei,ej,ek,sex,sey,sez,Imx,Jmy,Kmz;
1809: CellProperties cell_properties;
1810: PetscBool write_output = PETSC_FALSE;
1811: PetscErrorCode ierr;
1814: /* Generate the da for velocity and pressure */
1815: /* Num nodes in each direction is mx+1, my+1, mz+1 */
1816: u_dof = U_DOFS; /* Vx, Vy - velocities */
1817: p_dof = P_DOFS; /* p - pressure */
1818: dof = u_dof+p_dof;
1819: stencil_width = 1;
1820: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1821: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,PETSC_NULL,PETSC_NULL,PETSC_NULL,&da_Stokes);
1822: DMDASetFieldName(da_Stokes,0,"Vx");
1823: DMDASetFieldName(da_Stokes,1,"Vy");
1824: DMDASetFieldName(da_Stokes,2,"Vz");
1825: DMDASetFieldName(da_Stokes,3,"P");
1827: /* unit box [0,1] x [0,1] x [0,1] */
1828: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.0,1.0);
1830: /* local number of elements */
1831: DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,&mzl);
1833: /* create quadrature point info for PDE coefficients */
1834: CellPropertiesCreate(da_Stokes,&cell_properties);
1836: /* interpolate the coordinates to quadrature points */
1837: DMDAGetCoordinateDA(da_Stokes,&vel_cda);
1838: DMDAGetGhostedCoordinates(da_Stokes,&vel_coords);
1839: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1840: DMDAGetElementCorners(da_Stokes,&sex,&sey,&sez,&Imx,&Jmy,&Kmz);
1841: for (ek = sez; ek < sez+Kmz; ek++) {
1842: for (ej = sey; ej < sey+Jmy; ej++) {
1843: for (ei = sex; ei < sex+Imx; ei++) {
1844: /* get coords for the element */
1845: PetscInt ngp;
1846: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
1847: PetscScalar el_coords[NSD*NODES_PER_EL];
1848: GaussPointCoefficients *cell;
1850: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1851: GetElementCoords3D(_vel_coords,ei,ej,ek,el_coords);
1852: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1854: for (p = 0; p < GAUSS_POINTS; p++) {
1855: PetscScalar xi_p[NSD],Ni_p[NODES_PER_EL];
1856: PetscScalar gp_x,gp_y,gp_z;
1857: PetscInt n;
1859: xi_p[0] = gp_xi[p][0];
1860: xi_p[1] = gp_xi[p][1];
1861: xi_p[2] = gp_xi[p][2];
1862: ShapeFunctionQ13D_Evaluate(xi_p,Ni_p);
1864: gp_x = gp_y = gp_z = 0.0;
1865: for (n = 0; n < NODES_PER_EL; n++) {
1866: gp_x = gp_x+Ni_p[n]*el_coords[NSD*n ];
1867: gp_y = gp_y+Ni_p[n]*el_coords[NSD*n+1];
1868: gp_z = gp_z+Ni_p[n]*el_coords[NSD*n+2];
1869: }
1870: cell->gp_coords[NSD*p ] = gp_x;
1871: cell->gp_coords[NSD*p+1] = gp_y;
1872: cell->gp_coords[NSD*p+2] = gp_z;
1873: }
1874: }
1875: }
1876: }
1878: PetscOptionsGetInt(PETSC_NULL,"-model",&model_definition,PETSC_NULL);
1880: switch (model_definition) {
1881: case 0: /* isoviscous */
1882: for (ek = sez; ek < sez+Kmz; ek++) {
1883: for (ej = sey; ej < sey+Jmy; ej++) {
1884: for (ei = sex; ei < sex+Imx; ei++) {
1885: GaussPointCoefficients *cell;
1887: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1888: for (p = 0; p < GAUSS_POINTS; p++) {
1889: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p ]);
1890: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1891: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1893: cell->eta[p] = 1.0;
1895: cell->fx[p] = 0.0*coord_x;
1896: cell->fy[p] = -sin((double)2.2*M_PI*coord_y)*cos(1.0*M_PI*coord_x);
1897: cell->fz[p] = 0.0*coord_z;
1898: cell->hc[p] = 0.0;
1899: }
1900: }
1901: }
1902: }
1903: break;
1905: case 1: /* manufactured */
1906: for (ek = sez; ek < sez+Kmz; ek++) {
1907: for (ej = sey; ej < sey+Jmy; ej++) {
1908: for (ei = sex; ei < sex+Imx; ei++) {
1909: PetscReal eta,Fm[NSD],Fc,pos[NSD];
1910: GaussPointCoefficients *cell;
1912: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1913: for (p = 0; p < GAUSS_POINTS; p++) {
1914: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p ]);
1915: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1916: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1918: pos[0] = coord_x;
1919: pos[1] = coord_y;
1920: pos[2] = coord_z;
1922: evaluate_MS_FrankKamentski(pos,PETSC_NULL,PETSC_NULL,&eta,Fm,&Fc);
1923: cell->eta[p] = eta;
1925: cell->fx[p] = Fm[0];
1926: cell->fy[p] = Fm[1];
1927: cell->fz[p] = Fm[2];
1928: cell->hc[p] = Fc;
1929: }
1930: }
1931: }
1932: }
1933: break;
1935: case 2: /* solcx */
1936: for (ek = sez; ek < sez+Kmz; ek++) {
1937: for (ej = sey; ej < sey+Jmy; ej++) {
1938: for (ei = sex; ei < sex+Imx; ei++) {
1939: GaussPointCoefficients *cell;
1941: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1942: for (p = 0; p < GAUSS_POINTS; p++) {
1943: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p ]);
1944: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1945: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1947: cell->eta[p] = 1.0;
1949: cell->fx[p] = 0.0;
1950: cell->fy[p] = -sin((double)3*M_PI*coord_y)*cos(1.0*M_PI*coord_x);
1951: cell->fz[p] = 0.0*coord_z;
1952: cell->hc[p] = 0.0;
1953: }
1954: }
1955: }
1956: }
1957: break;
1959: case 3: /* sinker */
1960: for (ek = sez; ek < sez+Kmz; ek++) {
1961: for (ej = sey; ej < sey+Jmy; ej++) {
1962: for (ei = sex; ei < sex+Imx; ei++) {
1963: GaussPointCoefficients *cell;
1965: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1966: for (p = 0; p < GAUSS_POINTS; p++) {
1967: PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p ]);
1968: PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1969: PetscReal zp = PetscRealPart(cell->gp_coords[NSD*p+2]);
1971: cell->eta[p] = 1.0e-2;
1972: cell->fx[p] = 0.0;
1973: cell->fy[p] = 0.0;
1974: cell->fz[p] = 0.0;
1975: cell->hc[p] = 0.0;
1977: if ( (fabs(xp-0.5) < 0.2) &&
1978: (fabs(yp-0.5) < 0.2) &&
1979: (fabs(zp-0.5) < 0.2) ) {
1980: cell->eta[p] = 1.0;
1981: cell->fz[p] = 1.0;
1982: }
1984: }
1985: }
1986: }
1987: }
1988: break;
1990: default:
1991: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No default model is supported. Choose either -model {0,1,2,3}");
1992: break;
1993: }
1995: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1997: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1998: DMCreateMatrix(da_Stokes,MATAIJ,&A);
1999: DMCreateMatrix(da_Stokes,MATAIJ,&B);
2000: DMCreateGlobalVector(da_Stokes,&X);
2001: DMCreateGlobalVector(da_Stokes,&f);
2003: /* assemble A11 */
2004: MatZeroEntries(A);
2005: MatZeroEntries(B);
2006: VecZeroEntries(f);
2008: AssembleA_Stokes(A,da_Stokes,cell_properties);
2009: AssembleA_PCStokes(B,da_Stokes,cell_properties);
2010: /* build force vector */
2011: AssembleF_Stokes(f,da_Stokes,cell_properties);
2013: /* SOLVE */
2014: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
2015: KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
2016: KSPSetOperators(ksp_S,A,B,SAME_NONZERO_PATTERN);
2017: KSPSetFromOptions(ksp_S);
2019: {
2020: PC pc;
2021: const PetscInt ufields[] = {0,1,2},pfields[] = {3};
2022: KSPGetPC(ksp_S,&pc);
2023: PCFieldSplitSetBlockSize(pc,4);
2024: PCFieldSplitSetFields(pc,"u",3,ufields,ufields);
2025: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
2026: }
2028: {
2029: PC pc;
2030: PetscBool same = PETSC_FALSE;
2031: KSPGetPC(ksp_S,&pc);
2032: PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
2033: if (same) {
2034: PCMGSetupViaCoarsen(pc,da_Stokes);
2035: }
2036: }
2038: {
2039: PetscBool stokes_monitor = PETSC_FALSE;
2040: PetscOptionsGetBool(PETSC_NULL,"-stokes_ksp_monitor_blocks",&stokes_monitor,0);
2041: if (stokes_monitor) {
2042: KSPMonitorSet(ksp_S,KSPMonitorStokesBlocks,PETSC_NULL,PETSC_NULL);
2043: }
2044: }
2045: KSPSolve(ksp_S,f,X);
2047: PetscOptionsGetBool(PETSC_NULL,"-write_pvts",&write_output,PETSC_NULL);
2048: if (write_output) {DAView3DPVTS(da_Stokes,X,"up");}
2049: {
2050: PetscBool flg = PETSC_FALSE;
2051: char filename[PETSC_MAX_PATH_LEN];
2052: PetscOptionsGetString(PETSC_NULL,"-write_binary",filename,sizeof filename,&flg);
2053: if (flg) {
2054: PetscViewer viewer;
2055: //PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer);
2056: PetscViewerVTKOpen(PETSC_COMM_WORLD,"ex42.vts",FILE_MODE_WRITE,&viewer);
2057: VecView(X,viewer);
2058: PetscViewerDestroy(&viewer);
2059: }
2060: }
2061: KSPGetIterationNumber(ksp_S,&its);
2063: /* verify */
2064: if (model_definition == 1) {
2065: DM da_Stokes_analytic;
2066: Vec X_analytic;
2068: DMDACreateManufacturedSolution(mx,my,mz,&da_Stokes_analytic,&X_analytic);
2069: if (write_output) {
2070: DAView3DPVTS(da_Stokes_analytic,X_analytic,"ms");
2071: }
2072: DMDAIntegrateErrors3D(da_Stokes_analytic,X,X_analytic);
2073: if (write_output) {
2074: DAView3DPVTS(da_Stokes,X,"up2");
2075: }
2076: DMDestroy(&da_Stokes_analytic);
2077: VecDestroy(&X_analytic);
2078: }
2080: KSPDestroy(&ksp_S);
2081: VecDestroy(&X);
2082: VecDestroy(&f);
2083: MatDestroy(&A);
2084: MatDestroy(&B);
2086: CellPropertiesDestroy(&cell_properties);
2087: DMDestroy(&da_Stokes);
2088: return(0);
2089: }
2093: int main(int argc,char **args)
2094: {
2096: PetscInt mx,my,mz;
2098: PetscInitialize(&argc,&args,(char *)0,help);
2100: mx = my = mz = 10;
2101: PetscOptionsGetInt(PETSC_NULL,"-mx",&mx,PETSC_NULL);
2102: my = mx; mz = mx;
2103: PetscOptionsGetInt(PETSC_NULL,"-my",&my,PETSC_NULL);
2104: PetscOptionsGetInt(PETSC_NULL,"-mz",&mz,PETSC_NULL);
2106: solve_stokes_3d_coupled(mx,my,mz);
2108: PetscFinalize();
2109: return 0;
2110: }