Actual source code: ex42.c
petsc-3.4.5 2014-06-29
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) \n";
10: /* Contributed by Dave May */
12: #include petscksp.h
13: #include petscdmda.h
14: #include <petsctime.h>
16: #define PROFILE_TIMING
17: #define ASSEMBLE_LOWER_TRIANGULAR
19: #define NSD 3 /* number of spatial dimensions */
20: #define NODES_PER_EL 8 /* nodes per element */
21: #define U_DOFS 3 /* degrees of freedom per velocity node */
22: #define P_DOFS 1 /* degrees of freedom per pressure node */
23: #define GAUSS_POINTS 8
25: /* Gauss point based evaluation */
26: typedef struct {
27: PetscScalar gp_coords[NSD*GAUSS_POINTS];
28: PetscScalar eta[GAUSS_POINTS];
29: PetscScalar fx[GAUSS_POINTS];
30: PetscScalar fy[GAUSS_POINTS];
31: PetscScalar fz[GAUSS_POINTS];
32: PetscScalar hc[GAUSS_POINTS];
33: } GaussPointCoefficients;
35: typedef struct {
36: PetscScalar u_dof;
37: PetscScalar v_dof;
38: PetscScalar w_dof;
39: PetscScalar p_dof;
40: } StokesDOF;
42: typedef struct _p_CellProperties *CellProperties;
43: struct _p_CellProperties {
44: PetscInt ncells;
45: PetscInt mx,my,mz;
46: PetscInt sex,sey,sez;
47: GaussPointCoefficients *gpc;
48: };
50: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz);
52: /* elements */
55: PetscErrorCode CellPropertiesCreate(DM da_stokes,CellProperties *C)
56: {
58: CellProperties cells;
59: PetscInt mx,my,mz,sex,sey,sez;
62: PetscMalloc(sizeof(struct _p_CellProperties),&cells);
64: 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;
74: PetscMalloc(sizeof(GaussPointCoefficients)*mx*my*mz,&cells->gpc);
76: *C = cells;
77: return(0);
78: }
82: PetscErrorCode CellPropertiesDestroy(CellProperties *C)
83: {
85: CellProperties cells;
88: if (!C) return(0);
89: cells = *C;
90: PetscFree(cells->gpc);
91: PetscFree(cells);
92: *C = NULL;
93: return(0);
94: }
98: PetscErrorCode CellPropertiesGetCell(CellProperties C,PetscInt I,PetscInt J,PetscInt K,GaussPointCoefficients **G)
99: {
101: *G = &C->gpc[(I-C->sex) + (J-C->sey)*C->mx + (K-C->sez)*C->mx*C->my];
102: return(0);
103: }
105: /* FEM routines */
106: /*
107: Element: Local basis function ordering
108: 1-----2
109: | |
110: | |
111: 0-----3
112: */
113: static void ShapeFunctionQ13D_Evaluate(PetscScalar _xi[],PetscScalar Ni[])
114: {
115: PetscReal xi = PetscRealPart(_xi[0]);
116: PetscReal eta = PetscRealPart(_xi[1]);
117: PetscReal zeta = PetscRealPart(_xi[2]);
119: Ni[0] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 - zeta);
120: Ni[1] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 - zeta);
121: Ni[2] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 - zeta);
122: Ni[3] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 - zeta);
124: Ni[4] = 0.125 * (1.0 - xi) * (1.0 - eta) * (1.0 + zeta);
125: Ni[5] = 0.125 * (1.0 - xi) * (1.0 + eta) * (1.0 + zeta);
126: Ni[6] = 0.125 * (1.0 + xi) * (1.0 + eta) * (1.0 + zeta);
127: Ni[7] = 0.125 * (1.0 + xi) * (1.0 - eta) * (1.0 + zeta);
128: }
130: static void ShapeFunctionQ13D_Evaluate_dxi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])
131: {
132: PetscReal xi = PetscRealPart(_xi[0]);
133: PetscReal eta = PetscRealPart(_xi[1]);
134: PetscReal zeta = PetscRealPart(_xi[2]);
135: /* xi deriv */
136: GNi[0][0] = -0.125 * (1.0 - eta) * (1.0 - zeta);
137: GNi[0][1] = -0.125 * (1.0 + eta) * (1.0 - zeta);
138: GNi[0][2] = 0.125 * (1.0 + eta) * (1.0 - zeta);
139: GNi[0][3] = 0.125 * (1.0 - eta) * (1.0 - zeta);
141: GNi[0][4] = -0.125 * (1.0 - eta) * (1.0 + zeta);
142: GNi[0][5] = -0.125 * (1.0 + eta) * (1.0 + zeta);
143: GNi[0][6] = 0.125 * (1.0 + eta) * (1.0 + zeta);
144: GNi[0][7] = 0.125 * (1.0 - eta) * (1.0 + zeta);
145: /* eta deriv */
146: GNi[1][0] = -0.125 * (1.0 - xi) * (1.0 - zeta);
147: GNi[1][1] = 0.125 * (1.0 - xi) * (1.0 - zeta);
148: GNi[1][2] = 0.125 * (1.0 + xi) * (1.0 - zeta);
149: GNi[1][3] = -0.125 * (1.0 + xi) * (1.0 - zeta);
151: GNi[1][4] = -0.125 * (1.0 - xi) * (1.0 + zeta);
152: GNi[1][5] = 0.125 * (1.0 - xi) * (1.0 + zeta);
153: GNi[1][6] = 0.125 * (1.0 + xi) * (1.0 + zeta);
154: GNi[1][7] = -0.125 * (1.0 + xi) * (1.0 + zeta);
155: /* zeta deriv */
156: GNi[2][0] = -0.125 * (1.0 - xi) * (1.0 - eta);
157: GNi[2][1] = -0.125 * (1.0 - xi) * (1.0 + eta);
158: GNi[2][2] = -0.125 * (1.0 + xi) * (1.0 + eta);
159: GNi[2][3] = -0.125 * (1.0 + xi) * (1.0 - eta);
161: GNi[2][4] = 0.125 * (1.0 - xi) * (1.0 - eta);
162: GNi[2][5] = 0.125 * (1.0 - xi) * (1.0 + eta);
163: GNi[2][6] = 0.125 * (1.0 + xi) * (1.0 + eta);
164: GNi[2][7] = 0.125 * (1.0 + xi) * (1.0 - eta);
165: }
167: static void matrix_inverse_3x3(PetscScalar A[3][3],PetscScalar B[3][3])
168: {
169: PetscScalar t4, t6, t8, t10, t12, t14, t17;
171: t4 = A[2][0] * A[0][1];
172: t6 = A[2][0] * A[0][2];
173: t8 = A[1][0] * A[0][1];
174: t10 = A[1][0] * A[0][2];
175: t12 = A[0][0] * A[1][1];
176: t14 = A[0][0] * A[1][2];
177: 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]);
179: B[0][0] = (A[1][1] * A[2][2] - A[1][2] * A[2][1]) * t17;
180: B[0][1] = -(A[0][1] * A[2][2] - A[0][2] * A[2][1]) * t17;
181: B[0][2] = (A[0][1] * A[1][2] - A[0][2] * A[1][1]) * t17;
182: B[1][0] = -(-A[2][0] * A[1][2] + A[1][0] * A[2][2]) * t17;
183: B[1][1] = (-t6 + A[0][0] * A[2][2]) * t17;
184: B[1][2] = -(-t10 + t14) * t17;
185: B[2][0] = (-A[2][0] * A[1][1] + A[1][0] * A[2][1]) * t17;
186: B[2][1] = -(-t4 + A[0][0] * A[2][1]) * t17;
187: B[2][2] = (-t8 + t12) * t17;
188: }
190: static void ShapeFunctionQ13D_Evaluate_dx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
191: {
192: PetscScalar J00,J01,J02,J10,J11,J12,J20,J21,J22;
193: PetscInt n;
194: PetscScalar iJ[3][3],JJ[3][3];
196: J00 = J01 = J02 = 0.0;
197: J10 = J11 = J12 = 0.0;
198: J20 = J21 = J22 = 0.0;
199: for (n=0; n<NODES_PER_EL; n++) {
200: PetscScalar cx = coords[NSD*n + 0];
201: PetscScalar cy = coords[NSD*n + 1];
202: PetscScalar cz = coords[NSD*n + 2];
204: /* J_ij = d(x_j) / d(xi_i) */ /* J_ij = \sum _I GNi[j][I} * x_i */
205: J00 = J00 + GNi[0][n] * cx; /* J_xx */
206: J01 = J01 + GNi[0][n] * cy; /* J_xy = dx/deta */
207: J02 = J02 + GNi[0][n] * cz; /* J_xz = dx/dzeta */
209: J10 = J10 + GNi[1][n] * cx; /* J_yx = dy/dxi */
210: J11 = J11 + GNi[1][n] * cy; /* J_yy */
211: J12 = J12 + GNi[1][n] * cz; /* J_yz */
213: J20 = J20 + GNi[2][n] * cx; /* J_zx */
214: J21 = J21 + GNi[2][n] * cy; /* J_zy */
215: J22 = J22 + GNi[2][n] * cz; /* J_zz */
216: }
218: JJ[0][0] = J00; JJ[0][1] = J01; JJ[0][2] = J02;
219: JJ[1][0] = J10; JJ[1][1] = J11; JJ[1][2] = J12;
220: JJ[2][0] = J20; JJ[2][1] = J21; JJ[2][2] = J22;
222: matrix_inverse_3x3(JJ,iJ);
224: *det_J = J00*J11*J22 - J00*J12*J21 - J10*J01*J22 + J10*J02*J21 + J20*J01*J12 - J20*J02*J11;
226: for (n=0; n<NODES_PER_EL; n++) {
227: GNx[0][n] = GNi[0][n]*iJ[0][0] + GNi[1][n]*iJ[0][1] + GNi[2][n]*iJ[0][2];
228: GNx[1][n] = GNi[0][n]*iJ[1][0] + GNi[1][n]*iJ[1][1] + GNi[2][n]*iJ[1][2];
229: GNx[2][n] = GNi[0][n]*iJ[2][0] + GNi[1][n]*iJ[2][1] + GNi[2][n]*iJ[2][2];
230: }
231: }
233: static void ConstructGaussQuadrature3D(PetscInt *ngp,PetscScalar gp_xi[][NSD],PetscScalar gp_weight[])
234: {
235: *ngp = 8;
236: gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919; gp_xi[0][2] = -0.57735026919;
237: gp_xi[1][0] = -0.57735026919; gp_xi[1][1] = 0.57735026919; gp_xi[1][2] = -0.57735026919;
238: gp_xi[2][0] = 0.57735026919; gp_xi[2][1] = 0.57735026919; gp_xi[2][2] = -0.57735026919;
239: gp_xi[3][0] = 0.57735026919; gp_xi[3][1] = -0.57735026919; gp_xi[3][2] = -0.57735026919;
241: gp_xi[4][0] = -0.57735026919; gp_xi[4][1] = -0.57735026919; gp_xi[4][2] = 0.57735026919;
242: gp_xi[5][0] = -0.57735026919; gp_xi[5][1] = 0.57735026919; gp_xi[5][2] = 0.57735026919;
243: gp_xi[6][0] = 0.57735026919; gp_xi[6][1] = 0.57735026919; gp_xi[6][2] = 0.57735026919;
244: gp_xi[7][0] = 0.57735026919; gp_xi[7][1] = -0.57735026919; gp_xi[7][2] = 0.57735026919;
246: gp_weight[0] = 1.0;
247: gp_weight[1] = 1.0;
248: gp_weight[2] = 1.0;
249: gp_weight[3] = 1.0;
251: gp_weight[4] = 1.0;
252: gp_weight[5] = 1.0;
253: gp_weight[6] = 1.0;
254: gp_weight[7] = 1.0;
255: }
257: /* procs to the left claim the ghost node as their element */
260: static PetscErrorCode DMDAGetLocalElementSize(DM da,PetscInt *mxl,PetscInt *myl,PetscInt *mzl)
261: {
262: PetscInt m,n,p,M,N,P;
263: PetscInt sx,sy,sz;
267: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,0,0,0,0,0,0);
268: DMDAGetCorners(da,&sx,&sy,&sz,&m,&n,&p);
270: if (mxl) {
271: *mxl = m;
272: if ((sx+m) == M) *mxl = m-1; /* last proc */
273: }
274: if (myl) {
275: *myl = n;
276: if ((sy+n) == N) *myl = n-1; /* last proc */
277: }
278: if (mzl) {
279: *mzl = p;
280: if ((sz+p) == P) *mzl = p-1; /* last proc */
281: }
282: return(0);
283: }
287: static PetscErrorCode DMDAGetElementCorners(DM da,PetscInt *sx,PetscInt *sy,PetscInt *sz,PetscInt *mx,PetscInt *my,PetscInt *mz)
288: {
289: PetscInt si,sj,sk;
293: DMDAGetGhostCorners(da,&si,&sj,&sk,0,0,0);
295: if (sx) {
296: *sx = si;
297: if (si != 0) *sx = si+1;
298: }
299: if (sy) {
300: *sy = sj;
301: if (sj != 0) *sy = sj+1;
302: }
303: if (sz) {
304: *sz = sk;
305: if (sk != 0) *sz = sk+1;
306: }
307: DMDAGetLocalElementSize(da,mx,my,mz);
308: return(0);
309: }
311: /*
312: i,j are the element indices
313: The unknown is a vector quantity.
314: The s[].c is used to indicate the degree of freedom.
315: */
318: static PetscErrorCode DMDAGetElementEqnums3D_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j,PetscInt k)
319: {
320: PetscInt n;
323: /* velocity */
324: n = 0;
325: /* node 0 */
326: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++; /* Vx0 */
327: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++; /* Vy0 */
328: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++; /* Vz0 */
330: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
331: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
332: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
334: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 0; n++;
335: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 1; n++;
336: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k; s_u[n].c = 2; n++;
338: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 0; n++;
339: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 1; n++;
340: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k; s_u[n].c = 2; n++;
342: /* */
343: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++; /* Vx4 */
344: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++; /* Vy4 */
345: s_u[n].i = i; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++; /* Vz4 */
347: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
348: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
349: s_u[n].i = i; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
351: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 0; n++;
352: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 1; n++;
353: s_u[n].i = i+1; s_u[n].j = j+1; s_u[n].k = k+1; s_u[n].c = 2; n++;
355: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 0; n++;
356: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 1; n++;
357: s_u[n].i = i+1; s_u[n].j = j; s_u[n].k = k+1; s_u[n].c = 2; n++;
359: /* pressure */
360: n = 0;
362: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++; /* P0 */
363: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
364: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k; s_p[n].c = 3; n++;
365: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k; s_p[n].c = 3; n++;
367: s_p[n].i = i; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++; /* P0 */
368: s_p[n].i = i; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
369: s_p[n].i = i+1; s_p[n].j = j+1; s_p[n].k = k+1; s_p[n].c = 3; n++;
370: s_p[n].i = i+1; s_p[n].j = j; s_p[n].k = k+1; s_p[n].c = 3; n++;
371: return(0);
372: }
376: static PetscErrorCode GetElementCoords3D(DMDACoor3d ***coords,PetscInt i,PetscInt j,PetscInt k,PetscScalar el_coord[])
377: {
379: /* get coords for the element */
380: el_coord[0] = coords[k][j][i].x;
381: el_coord[1] = coords[k][j][i].y;
382: el_coord[2] = coords[k][j][i].z;
384: el_coord[3] = coords[k][j+1][i].x;
385: el_coord[4] = coords[k][j+1][i].y;
386: el_coord[5] = coords[k][j+1][i].z;
388: el_coord[6] = coords[k][j+1][i+1].x;
389: el_coord[7] = coords[k][j+1][i+1].y;
390: el_coord[8] = coords[k][j+1][i+1].z;
392: el_coord[9] = coords[k][j][i+1].x;
393: el_coord[10] = coords[k][j][i+1].y;
394: el_coord[11] = coords[k][j][i+1].z;
396: el_coord[12] = coords[k+1][j][i].x;
397: el_coord[13] = coords[k+1][j][i].y;
398: el_coord[14] = coords[k+1][j][i].z;
400: el_coord[15] = coords[k+1][j+1][i].x;
401: el_coord[16] = coords[k+1][j+1][i].y;
402: el_coord[17] = coords[k+1][j+1][i].z;
404: el_coord[18] = coords[k+1][j+1][i+1].x;
405: el_coord[19] = coords[k+1][j+1][i+1].y;
406: el_coord[20] = coords[k+1][j+1][i+1].z;
408: el_coord[21] = coords[k+1][j][i+1].x;
409: el_coord[22] = coords[k+1][j][i+1].y;
410: el_coord[23] = coords[k+1][j][i+1].z;
411: return(0);
412: }
416: static PetscErrorCode StokesDAGetNodalFields3D(StokesDOF ***field,PetscInt i,PetscInt j,PetscInt k,StokesDOF nodal_fields[])
417: {
419: /* get the nodal fields for u */
420: nodal_fields[0].u_dof = field[k][j][i].u_dof;
421: nodal_fields[0].v_dof = field[k][j][i].v_dof;
422: nodal_fields[0].w_dof = field[k][j][i].w_dof;
424: nodal_fields[1].u_dof = field[k][j+1][i].u_dof;
425: nodal_fields[1].v_dof = field[k][j+1][i].v_dof;
426: nodal_fields[1].w_dof = field[k][j+1][i].w_dof;
428: nodal_fields[2].u_dof = field[k][j+1][i+1].u_dof;
429: nodal_fields[2].v_dof = field[k][j+1][i+1].v_dof;
430: nodal_fields[2].w_dof = field[k][j+1][i+1].w_dof;
432: nodal_fields[3].u_dof = field[k][j][i+1].u_dof;
433: nodal_fields[3].v_dof = field[k][j][i+1].v_dof;
434: nodal_fields[3].w_dof = field[k][j][i+1].w_dof;
436: nodal_fields[4].u_dof = field[k+1][j][i].u_dof;
437: nodal_fields[4].v_dof = field[k+1][j][i].v_dof;
438: nodal_fields[4].w_dof = field[k+1][j][i].w_dof;
440: nodal_fields[5].u_dof = field[k+1][j+1][i].u_dof;
441: nodal_fields[5].v_dof = field[k+1][j+1][i].v_dof;
442: nodal_fields[5].w_dof = field[k+1][j+1][i].w_dof;
444: nodal_fields[6].u_dof = field[k+1][j+1][i+1].u_dof;
445: nodal_fields[6].v_dof = field[k+1][j+1][i+1].v_dof;
446: nodal_fields[6].w_dof = field[k+1][j+1][i+1].w_dof;
448: nodal_fields[7].u_dof = field[k+1][j][i+1].u_dof;
449: nodal_fields[7].v_dof = field[k+1][j][i+1].v_dof;
450: nodal_fields[7].w_dof = field[k+1][j][i+1].w_dof;
452: /* get the nodal fields for p */
453: nodal_fields[0].p_dof = field[k][j][i].p_dof;
454: nodal_fields[1].p_dof = field[k][j+1][i].p_dof;
455: nodal_fields[2].p_dof = field[k][j+1][i+1].p_dof;
456: nodal_fields[3].p_dof = field[k][j][i+1].p_dof;
458: nodal_fields[4].p_dof = field[k+1][j][i].p_dof;
459: nodal_fields[5].p_dof = field[k+1][j+1][i].p_dof;
460: nodal_fields[6].p_dof = field[k+1][j+1][i+1].p_dof;
461: nodal_fields[7].p_dof = field[k+1][j][i+1].p_dof;
462: return(0);
463: }
465: 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)
466: {
467: PetscInt ij;
468: PETSC_UNUSED PetscInt r,c,nr,nc;
470: nr = w_NPE*w_dof;
471: nc = u_NPE*u_dof;
473: r = w_dof*wi+wd;
474: c = u_dof*ui+ud;
476: ij = r*nc+c;
478: return ij;
479: }
483: static PetscErrorCode DMDASetValuesLocalStencil3D_ADD_VALUES(StokesDOF ***fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])
484: {
485: PetscInt n,I,J,K;
488: for (n = 0; n<NODES_PER_EL; n++) {
489: I = u_eqn[NSD*n].i;
490: J = u_eqn[NSD*n].j;
491: K = u_eqn[NSD*n].k;
493: fields_F[K][J][I].u_dof = fields_F[K][J][I].u_dof+Fe_u[NSD*n];
495: I = u_eqn[NSD*n+1].i;
496: J = u_eqn[NSD*n+1].j;
497: K = u_eqn[NSD*n+1].k;
499: fields_F[K][J][I].v_dof = fields_F[K][J][I].v_dof+Fe_u[NSD*n+1];
501: I = u_eqn[NSD*n+2].i;
502: J = u_eqn[NSD*n+2].j;
503: K = u_eqn[NSD*n+2].k;
504: fields_F[K][J][I].w_dof = fields_F[K][J][I].w_dof+Fe_u[NSD*n+2];
506: I = p_eqn[n].i;
507: J = p_eqn[n].j;
508: K = p_eqn[n].k;
510: fields_F[K][J][I].p_dof = fields_F[K][J][I].p_dof+Fe_p[n];
512: }
513: return(0);
514: }
516: static void FormStressOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
517: {
518: PetscInt ngp;
519: PetscScalar gp_xi[GAUSS_POINTS][NSD];
520: PetscScalar gp_weight[GAUSS_POINTS];
521: PetscInt p,i,j,k;
522: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
523: PetscScalar J_p,tildeD[6];
524: PetscScalar B[6][U_DOFS*NODES_PER_EL];
525: const PetscInt nvdof = U_DOFS*NODES_PER_EL;
527: /* define quadrature rule */
528: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
530: /* evaluate integral */
531: for (p = 0; p < ngp; p++) {
532: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
533: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
535: for (i = 0; i < NODES_PER_EL; i++) {
536: PetscScalar d_dx_i = GNx_p[0][i];
537: PetscScalar d_dy_i = GNx_p[1][i];
538: PetscScalar d_dz_i = GNx_p[2][i];
540: B[0][3*i] = d_dx_i; B[0][3*i+1] = 0.0; B[0][3*i+2] = 0.0;
541: B[1][3*i] = 0.0; B[1][3*i+1] = d_dy_i; B[1][3*i+2] = 0.0;
542: B[2][3*i] = 0.0; B[2][3*i+1] = 0.0; B[2][3*i+2] = d_dz_i;
544: B[3][3*i] = d_dy_i; B[3][3*i+1] = d_dx_i; B[3][3*i+2] = 0.0; /* e_xy */
545: B[4][3*i] = d_dz_i; B[4][3*i+1] = 0.0; B[4][3*i+2] = d_dx_i; /* e_xz */
546: B[5][3*i] = 0.0; B[5][3*i+1] = d_dz_i; B[5][3*i+2] = d_dy_i; /* e_yz */
547: }
550: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
551: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
552: tildeD[2] = 2.0*gp_weight[p]*J_p*eta[p];
554: tildeD[3] = gp_weight[p]*J_p*eta[p];
555: tildeD[4] = gp_weight[p]*J_p*eta[p];
556: tildeD[5] = gp_weight[p]*J_p*eta[p];
558: /* form Bt tildeD B */
559: /*
560: Ke_ij = Bt_ik . D_kl . B_lj
561: = B_ki . D_kl . B_lj
562: = B_ki . D_kk . B_kj
563: */
564: for (i = 0; i < nvdof; i++) {
565: for (j = i; j < nvdof; j++) {
566: for (k = 0; k < 6; k++) {
567: Ke[i*nvdof+j] += B[k][i]*tildeD[k]*B[k][j];
568: }
569: }
570: }
572: }
573: /* fill lower triangular part */
574: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
575: for (i = 0; i < nvdof; i++) {
576: for (j = i; j < nvdof; j++) {
577: Ke[j*nvdof+i] = Ke[i*nvdof+j];
578: }
579: }
580: #endif
581: }
583: static void FormGradientOperatorQ13D(PetscScalar Ke[],PetscScalar coords[])
584: {
585: PetscInt ngp;
586: PetscScalar gp_xi[GAUSS_POINTS][NSD];
587: PetscScalar gp_weight[GAUSS_POINTS];
588: PetscInt p,i,j,di;
589: PetscScalar Ni_p[NODES_PER_EL];
590: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
591: PetscScalar J_p,fac;
593: /* define quadrature rule */
594: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
596: /* evaluate integral */
597: for (p = 0; p < ngp; p++) {
598: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
599: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
600: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
601: fac = gp_weight[p]*J_p;
603: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
604: for (di = 0; di < NSD; di++) { /* u dofs */
605: for (j = 0; j < NODES_PER_EL; j++) { /* p nodes, p dofs = 1 (ie no loop) */
606: PetscInt IJ;
607: IJ = ASS_MAP_wIwDI_uJuDJ(i,di,NODES_PER_EL,3,j,0,NODES_PER_EL,1);
609: Ke[IJ] = Ke[IJ]-GNx_p[di][i]*Ni_p[j]*fac;
610: }
611: }
612: }
613: }
614: }
616: static void FormDivergenceOperatorQ13D(PetscScalar De[],PetscScalar coords[])
617: {
618: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
619: PetscInt i,j;
620: PetscInt nr_g,nc_g;
622: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
623: FormGradientOperatorQ13D(Ge,coords);
625: nr_g = U_DOFS*NODES_PER_EL;
626: nc_g = P_DOFS*NODES_PER_EL;
628: for (i = 0; i < nr_g; i++) {
629: for (j = 0; j < nc_g; j++) {
630: De[nr_g*j+i] = Ge[nc_g*i+j];
631: }
632: }
633: }
635: static void FormStabilisationOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
636: {
637: PetscInt ngp;
638: PetscScalar gp_xi[GAUSS_POINTS][NSD];
639: PetscScalar gp_weight[GAUSS_POINTS];
640: PetscInt p,i,j;
641: PetscScalar Ni_p[NODES_PER_EL];
642: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
643: PetscScalar J_p,fac,eta_avg;
645: /* define quadrature rule */
646: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
648: /* evaluate integral */
649: for (p = 0; p < ngp; p++) {
650: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
651: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
652: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
653: fac = gp_weight[p]*J_p;
654: /*
655: for (i = 0; i < NODES_PER_EL; i++) {
656: for (j = i; j < NODES_PER_EL; j++) {
657: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
658: }
659: }
660: */
662: for (i = 0; i < NODES_PER_EL; i++) {
663: for (j = 0; j < NODES_PER_EL; j++) {
664: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.015625);
665: }
666: }
667: }
669: /* scale */
670: eta_avg = 0.0;
671: for (p = 0; p < ngp; p++) eta_avg += eta[p];
672: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
673: fac = 1.0/eta_avg;
675: /*
676: for (i = 0; i < NODES_PER_EL; i++) {
677: for (j = i; j < NODES_PER_EL; j++) {
678: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
679: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
680: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
681: #endif
682: }
683: }
684: */
686: for (i = 0; i < NODES_PER_EL; i++) {
687: for (j = 0; j < NODES_PER_EL; j++) {
688: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
689: }
690: }
691: }
693: static void FormScaledMassMatrixOperatorQ13D(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
694: {
695: PetscInt ngp;
696: PetscScalar gp_xi[GAUSS_POINTS][NSD];
697: PetscScalar gp_weight[GAUSS_POINTS];
698: PetscInt p,i,j;
699: PetscScalar Ni_p[NODES_PER_EL];
700: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
701: PetscScalar J_p,fac,eta_avg;
703: /* define quadrature rule */
704: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
706: /* evaluate integral */
707: for (p = 0; p < ngp; p++) {
708: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
709: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
710: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
711: fac = gp_weight[p]*J_p;
713: /*
714: for (i = 0; i < NODES_PER_EL; i++) {
715: for (j = i; j < NODES_PER_EL; j++) {
716: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
717: }
718: }
719: */
721: for (i = 0; i < NODES_PER_EL; i++) {
722: for (j = 0; j < NODES_PER_EL; j++) {
723: Ke[NODES_PER_EL*i+j] = Ke[NODES_PER_EL*i+j]-fac*Ni_p[i]*Ni_p[j];
724: }
725: }
726: }
728: /* scale */
729: eta_avg = 0.0;
730: for (p = 0; p < ngp; p++) eta_avg += eta[p];
731: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
732: fac = 1.0/eta_avg;
733: /*
734: for (i = 0; i < NODES_PER_EL; i++) {
735: for (j = i; j < NODES_PER_EL; j++) {
736: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
737: #if defined(ASSEMBLE_LOWER_TRIANGULAR)
738: Ke[NODES_PER_EL*j+i] = Ke[NODES_PER_EL*i+j];
739: #endif
740: }
741: }
742: */
744: for (i = 0; i < NODES_PER_EL; i++) {
745: for (j = 0; j < NODES_PER_EL; j++) {
746: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
747: }
748: }
749: }
751: static void FormMomentumRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[],PetscScalar fz[])
752: {
753: PetscInt ngp;
754: PetscScalar gp_xi[GAUSS_POINTS][NSD];
755: PetscScalar gp_weight[GAUSS_POINTS];
756: PetscInt p,i;
757: PetscScalar Ni_p[NODES_PER_EL];
758: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
759: PetscScalar J_p,fac;
761: /* define quadrature rule */
762: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
764: /* evaluate integral */
765: for (p = 0; p < ngp; p++) {
766: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
767: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
768: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
769: fac = gp_weight[p]*J_p;
771: for (i = 0; i < NODES_PER_EL; i++) {
772: Fe[NSD*i] -= fac*Ni_p[i]*fx[p];
773: Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
774: Fe[NSD*i+2] -= fac*Ni_p[i]*fz[p];
775: }
776: }
777: }
779: static void FormContinuityRhsQ13D(PetscScalar Fe[],PetscScalar coords[],PetscScalar hc[])
780: {
781: PetscInt ngp;
782: PetscScalar gp_xi[GAUSS_POINTS][NSD];
783: PetscScalar gp_weight[GAUSS_POINTS];
784: PetscInt p,i;
785: PetscScalar Ni_p[NODES_PER_EL];
786: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
787: PetscScalar J_p,fac;
789: /* define quadrature rule */
790: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
792: /* evaluate integral */
793: for (p = 0; p < ngp; p++) {
794: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
795: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
796: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,coords,&J_p);
797: fac = gp_weight[p]*J_p;
799: for (i = 0; i < NODES_PER_EL; i++) Fe[i] -= fac*Ni_p[i]*hc[p];
800: }
801: }
803: #define _ZERO_ROWCOL_i(A,i) { \
804: PetscInt KK; \
805: PetscScalar tmp = A[24*(i)+(i)]; \
806: for (KK=0;KK<24;KK++) A[24*(i)+KK]=0.0; \
807: for (KK=0;KK<24;KK++) A[24*KK+(i)]=0.0; \
808: A[24*(i)+(i)] = tmp;} \
810: #define _ZERO_ROW_i(A,i) { \
811: PetscInt KK; \
812: for (KK=0;KK<8;KK++) A[8*(i)+KK]=0.0;}
814: #define _ZERO_COL_i(A,i) { \
815: PetscInt KK; \
816: for (KK=0;KK<8;KK++) A[24*KK+(i)]=0.0;}
820: static PetscErrorCode AssembleA_Stokes(Mat A,DM stokes_da,CellProperties cell_properties)
821: {
822: DM cda;
823: Vec coords;
824: DMDACoor3d ***_coords;
825: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
826: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
827: PetscInt sex,sey,sez,mx,my,mz;
828: PetscInt ei,ej,ek;
829: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
830: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
831: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
832: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
833: PetscScalar el_coords[NODES_PER_EL*NSD];
834: GaussPointCoefficients *props;
835: PetscScalar *prop_eta;
836: PetscInt n,M,N,P;
837: PetscLogDouble t0,t1;
838: PetscErrorCode ierr;
841: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
842: /* setup for coords */
843: DMGetCoordinateDM(stokes_da,&cda);
844: DMGetCoordinatesLocal(stokes_da,&coords);
845: DMDAVecGetArray(cda,coords,&_coords);
847: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
848: PetscTime(&t0);
849: for (ek = sez; ek < sez+mz; ek++) {
850: for (ej = sey; ej < sey+my; ej++) {
851: for (ei = sex; ei < sex+mx; ei++) {
852: /* get coords for the element */
853: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
854: /* get cell properties */
855: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
856: /* get coefficients for the element */
857: prop_eta = props->eta;
859: /* initialise element stiffness matrix */
860: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
861: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
862: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
863: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
865: /* form element stiffness matrix */
866: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
867: FormGradientOperatorQ13D(Ge,el_coords);
868: /*#if defined(ASSEMBLE_LOWER_TRIANGULAR)*/
869: FormDivergenceOperatorQ13D(De,el_coords);
870: /*#endif*/
871: FormStabilisationOperatorQ13D(Ce,el_coords,prop_eta);
873: /* insert element matrix into global matrix */
874: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
876: for (n=0; n<NODES_PER_EL; n++) {
877: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
878: _ZERO_ROWCOL_i(Ae,3*n);
879: _ZERO_ROW_i (Ge,3*n);
880: _ZERO_COL_i (De,3*n);
881: }
883: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
884: _ZERO_ROWCOL_i(Ae,3*n+1);
885: _ZERO_ROW_i (Ge,3*n+1);
886: _ZERO_COL_i (De,3*n+1);
887: }
889: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
890: _ZERO_ROWCOL_i(Ae,3*n+2);
891: _ZERO_ROW_i (Ge,3*n+2);
892: _ZERO_COL_i (De,3*n+2);
893: }
894: }
895: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
896: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
897: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
898: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
899: }
900: }
901: }
902: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
903: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
905: DMDAVecRestoreArray(cda,coords,&_coords);
907: PetscTime(&t1);
908: return(0);
909: }
913: static PetscErrorCode AssembleA_PCStokes(Mat A,DM stokes_da,CellProperties cell_properties)
914: {
915: DM cda;
916: Vec coords;
917: DMDACoor3d ***_coords;
918: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
919: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
920: PetscInt sex,sey,sez,mx,my,mz;
921: PetscInt ei,ej,ek;
922: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
923: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
924: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
925: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
926: PetscScalar el_coords[NODES_PER_EL*NSD];
927: GaussPointCoefficients *props;
928: PetscScalar *prop_eta;
929: PetscInt n,M,N,P;
930: PetscErrorCode ierr;
933: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
934: /* setup for coords */
935: DMGetCoordinateDM(stokes_da,&cda);
936: DMGetCoordinatesLocal(stokes_da,&coords);
937: DMDAVecGetArray(cda,coords,&_coords);
939: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
940: for (ek = sez; ek < sez+mz; ek++) {
941: for (ej = sey; ej < sey+my; ej++) {
942: for (ei = sex; ei < sex+mx; ei++) {
943: /* get coords for the element */
944: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
945: /* get cell properties */
946: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
947: /* get coefficients for the element */
948: prop_eta = props->eta;
950: /* initialise element stiffness matrix */
951: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
952: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
953: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
954: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
956: /* form element stiffness matrix */
957: FormStressOperatorQ13D(Ae,el_coords,prop_eta);
958: FormGradientOperatorQ13D(Ge,el_coords);
959: /* FormDivergenceOperatorQ13D(De,el_coords); */
960: FormScaledMassMatrixOperatorQ13D(Ce,el_coords,prop_eta);
962: /* insert element matrix into global matrix */
963: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
965: for (n=0; n<NODES_PER_EL; n++) {
966: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) {
967: _ZERO_ROWCOL_i(Ae,3*n);
968: _ZERO_ROW_i (Ge,3*n);
969: _ZERO_COL_i (De,3*n);
970: }
972: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) {
973: _ZERO_ROWCOL_i(Ae,3*n+1);
974: _ZERO_ROW_i (Ge,3*n+1);
975: _ZERO_COL_i (De,3*n+1);
976: }
978: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) {
979: _ZERO_ROWCOL_i(Ae,3*n+2);
980: _ZERO_ROW_i (Ge,3*n+2);
981: _ZERO_COL_i (De,3*n+2);
982: }
983: }
984: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
985: MatSetValuesStencil(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
986: /*MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);*/
987: MatSetValuesStencil(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
988: }
989: }
990: }
991: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
992: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
994: DMDAVecRestoreArray(cda,coords,&_coords);
995: return(0);
996: }
1000: static PetscErrorCode AssembleF_Stokes(Vec F,DM stokes_da,CellProperties cell_properties)
1001: {
1002: DM cda;
1003: Vec coords;
1004: DMDACoor3d ***_coords;
1005: MatStencil u_eqn[NODES_PER_EL*U_DOFS];
1006: MatStencil p_eqn[NODES_PER_EL*P_DOFS];
1007: PetscInt sex,sey,sez,mx,my,mz;
1008: PetscInt ei,ej,ek;
1009: PetscScalar Fe[NODES_PER_EL*U_DOFS];
1010: PetscScalar He[NODES_PER_EL*P_DOFS];
1011: PetscScalar el_coords[NODES_PER_EL*NSD];
1012: GaussPointCoefficients *props;
1013: PetscScalar *prop_fx,*prop_fy,*prop_fz,*prop_hc;
1014: Vec local_F;
1015: StokesDOF ***ff;
1016: PetscInt n,M,N,P;
1017: PetscErrorCode ierr;
1020: DMDAGetInfo(stokes_da,0,&M,&N,&P,0,0,0, 0,0,0,0,0,0);
1021: /* setup for coords */
1022: DMGetCoordinateDM(stokes_da,&cda);
1023: DMGetCoordinatesLocal(stokes_da,&coords);
1024: DMDAVecGetArray(cda,coords,&_coords);
1026: /* get acces to the vector */
1027: DMGetLocalVector(stokes_da,&local_F);
1028: VecZeroEntries(local_F);
1029: DMDAVecGetArray(stokes_da,local_F,&ff);
1031: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1032: for (ek = sez; ek < sez+mz; ek++) {
1033: for (ej = sey; ej < sey+my; ej++) {
1034: for (ei = sex; ei < sex+mx; ei++) {
1035: /* get coords for the element */
1036: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1037: /* get cell properties */
1038: CellPropertiesGetCell(cell_properties,ei,ej,ek,&props);
1039: /* get coefficients for the element */
1040: prop_fx = props->fx;
1041: prop_fy = props->fy;
1042: prop_fz = props->fz;
1043: prop_hc = props->hc;
1045: /* initialise element stiffness matrix */
1046: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
1047: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
1049: /* form element stiffness matrix */
1050: FormMomentumRhsQ13D(Fe,el_coords,prop_fx,prop_fy,prop_fz);
1051: FormContinuityRhsQ13D(He,el_coords,prop_hc);
1053: /* insert element matrix into global matrix */
1054: DMDAGetElementEqnums3D_up(u_eqn,p_eqn,ei,ej,ek);
1056: for (n=0; n<NODES_PER_EL; n++) {
1057: if ((u_eqn[3*n].i == 0) || (u_eqn[3*n].i == M-1)) Fe[3*n] = 0.0;
1059: if ((u_eqn[3*n+1].j == 0) || (u_eqn[3*n+1].j == N-1)) Fe[3*n+1] = 0.0;
1061: if ((u_eqn[3*n+2].k == 0) || (u_eqn[3*n+2].k == P-1)) Fe[3*n+2] = 0.0;
1062: }
1064: DMDASetValuesLocalStencil3D_ADD_VALUES(ff,u_eqn,p_eqn,Fe,He);
1065: }
1066: }
1067: }
1068: DMDAVecRestoreArray(stokes_da,local_F,&ff);
1069: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
1070: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
1071: DMRestoreLocalVector(stokes_da,&local_F);
1073: DMDAVecRestoreArray(cda,coords,&_coords);
1074: return(0);
1075: }
1077: static void evaluate_MS_FrankKamentski_constants(PetscReal *theta,PetscReal *MX,PetscReal *MY,PetscReal *MZ)
1078: {
1079: *theta = 0.0;
1080: *MX = 2.0 * M_PI;
1081: *MY = 2.0 * M_PI;
1082: *MZ = 2.0 * M_PI;
1083: }
1084: static void evaluate_MS_FrankKamentski(PetscReal pos[],PetscReal v[],PetscReal *p,PetscReal *eta,PetscReal Fm[],PetscReal *Fc)
1085: {
1086: PetscReal x,y,z;
1087: PetscReal theta,MX,MY,MZ;
1089: evaluate_MS_FrankKamentski_constants(&theta,&MX,&MY,&MZ);
1090: x = pos[0];
1091: y = pos[1];
1092: z = pos[2];
1093: if (v) {
1094: /*
1095: v[0] = pow(z,3)*exp(y)*sin(M_PI*x);
1096: v[1] = z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y);
1097: v[2] = -(pow(x,3) + pow(y,3))*sin(2.0*M_PI*z);
1098: */
1099: /*
1100: v[0] = sin(M_PI*x);
1101: v[1] = sin(M_PI*y);
1102: v[2] = sin(M_PI*z);
1103: */
1104: v[0] = PetscPowRealInt(z,3)*exp(y)*sin(M_PI*x);
1105: v[1] = z*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y);
1106: v[2] = PetscPowRealInt(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*PetscPowRealInt(z,4)*cos(M_PI*x)*exp(y)/4;
1107: }
1108: if (p) *p = PetscPowRealInt(x,2) + PetscPowRealInt(y,2) + PetscPowRealInt(z,2);
1109: if (eta) {
1110: /**eta = exp(-theta*(1.0 - y - 0.1*cos(MX*x)*cos(MZ*z)*sin(MY*y)));*/
1111: *eta = 1.0;
1112: }
1113: if (Fm) {
1114: /*
1115: 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))) ;
1116: 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)));
1117: 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))) ;
1118: */
1119: /*
1120: Fm[0]=-2*x - 2.0*pow(M_PI,2)*sin(M_PI*x);
1121: Fm[1]=-2*y - 2.0*pow(M_PI,2)*sin(M_PI*y);
1122: Fm[2]=-2*z - 2.0*pow(M_PI,2)*sin(M_PI*z);
1123: */
1124: /*
1125: 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) ;
1126: 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);
1127: 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) ;
1128: */
1129: Fm[0] = -2*x + 2*z*(PetscPowRealInt(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)) + PetscPowRealInt(z,3)*exp(y)*sin(M_PI*x) + 6.0*z*exp(y)*sin(M_PI*x) - 1.0*PetscPowRealInt(M_PI,2)*PetscPowRealInt(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*PetscPowRealInt(M_PI,2)*cos(M_PI*y)*exp(-y)*sin(2.0*M_PI*x);
1130: Fm[1] = -2*y + 2*z*(-cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + PetscPowRealInt(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*PetscPowRealInt(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);
1131: Fm[2] = -2*z + PetscPowRealInt(z,2)*(-2.0*PetscPowRealInt(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y) + 2.0*PetscPowRealInt(M_PI,3)*cos(M_PI*y)*cos(2.0*M_PI*x)*exp(-y)) + PetscPowRealInt(z,2)*(cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 - 3*PetscPowRealInt(M_PI,2)*cos(2.0*M_PI*x)*exp(-y)*sin(M_PI*y)/2 + PetscPowRealInt(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*PetscPowRealInt(M_PI,3)*PetscPowRealInt(z,4)*cos(M_PI*x)*exp(y) - 0.25*M_PI*PetscPowRealInt(z,4)*cos(M_PI*x)*exp(y) - 3.0*M_PI*PetscPowRealInt(z,2)*cos(M_PI*x)*exp(y) - 1.0*M_PI*cos(M_PI *y)*cos(2.0*M_PI*x)*exp(-y);
1132: }
1133: if (Fc) {
1134: /**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) ;*/
1135: /**Fc = M_PI*cos(M_PI*x) + M_PI*cos(M_PI*y) + M_PI*cos(M_PI*z);*/
1136: /**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);*/
1137: *Fc = 0.0;
1138: }
1139: }
1143: static PetscErrorCode DMDACreateManufacturedSolution(PetscInt mx,PetscInt my,PetscInt mz,DM *_da,Vec *_X)
1144: {
1145: DM da,cda;
1146: Vec X,local_X;
1147: StokesDOF ***_stokes;
1148: Vec coords;
1149: DMDACoor3d ***_coords;
1150: PetscInt si,sj,sk,ei,ej,ek,i,j,k;
1154: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1155: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,4,1,NULL,NULL,NULL,&da);
1156: DMDASetFieldName(da,0,"anlytic_Vx");
1157: DMDASetFieldName(da,1,"anlytic_Vy");
1158: DMDASetFieldName(da,2,"anlytic_Vz");
1159: DMDASetFieldName(da,3,"analytic_P");
1161: DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,1.0);
1163: DMGetCoordinatesLocal(da,&coords);
1164: DMGetCoordinateDM(da,&cda);
1165: DMDAVecGetArray(cda,coords,&_coords);
1167: DMCreateGlobalVector(da,&X);
1168: DMCreateLocalVector(da,&local_X);
1169: DMDAVecGetArray(da,local_X,&_stokes);
1171: DMDAGetGhostCorners(da,&si,&sj,&sk,&ei,&ej,&ek);
1172: for (k = sk; k < sk+ek; k++) {
1173: for (j = sj; j < sj+ej; j++) {
1174: for (i = si; i < si+ei; i++) {
1175: PetscReal pos[NSD],pressure,vel[NSD];
1177: pos[0] = PetscRealPart(_coords[k][j][i].x);
1178: pos[1] = PetscRealPart(_coords[k][j][i].y);
1179: pos[2] = PetscRealPart(_coords[k][j][i].z);
1181: evaluate_MS_FrankKamentski(pos,vel,&pressure,NULL,NULL,NULL);
1183: _stokes[k][j][i].u_dof = vel[0];
1184: _stokes[k][j][i].v_dof = vel[1];
1185: _stokes[k][j][i].w_dof = vel[2];
1186: _stokes[k][j][i].p_dof = pressure;
1187: }
1188: }
1189: }
1190: DMDAVecRestoreArray(da,local_X,&_stokes);
1191: DMDAVecRestoreArray(cda,coords,&_coords);
1193: DMLocalToGlobalBegin(da,local_X,INSERT_VALUES,X);
1194: DMLocalToGlobalEnd(da,local_X,INSERT_VALUES,X);
1196: VecDestroy(&local_X);
1198: *_da = da;
1199: *_X = X;
1200: return(0);
1201: }
1205: static PetscErrorCode DMDAIntegrateErrors3D(DM stokes_da,Vec X,Vec X_analytic)
1206: {
1207: DM cda;
1208: Vec coords,X_analytic_local,X_local;
1209: DMDACoor3d ***_coords;
1210: PetscInt sex,sey,sez,mx,my,mz;
1211: PetscInt ei,ej,ek;
1212: PetscScalar el_coords[NODES_PER_EL*NSD];
1213: StokesDOF ***stokes_analytic,***stokes;
1214: StokesDOF stokes_analytic_e[NODES_PER_EL],stokes_e[NODES_PER_EL];
1216: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
1217: PetscScalar Ni_p[NODES_PER_EL];
1218: PetscInt ngp;
1219: PetscScalar gp_xi[GAUSS_POINTS][NSD];
1220: PetscScalar gp_weight[GAUSS_POINTS];
1221: PetscInt p,i;
1222: PetscScalar J_p,fac;
1223: PetscScalar h,p_e_L2,u_e_L2,u_e_H1,p_L2,u_L2,u_H1,tp_L2,tu_L2,tu_H1;
1224: PetscScalar tint_p_ms,tint_p,int_p_ms,int_p;
1225: PetscInt M;
1226: PetscReal xymin[NSD],xymax[NSD];
1230: /* define quadrature rule */
1231: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1233: /* setup for coords */
1234: DMGetCoordinateDM(stokes_da,&cda);
1235: DMGetCoordinatesLocal(stokes_da,&coords);
1236: DMDAVecGetArray(cda,coords,&_coords);
1238: /* setup for analytic */
1239: DMCreateLocalVector(stokes_da,&X_analytic_local);
1240: DMGlobalToLocalBegin(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1241: DMGlobalToLocalEnd(stokes_da,X_analytic,INSERT_VALUES,X_analytic_local);
1242: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1244: /* setup for solution */
1245: DMCreateLocalVector(stokes_da,&X_local);
1246: DMGlobalToLocalBegin(stokes_da,X,INSERT_VALUES,X_local);
1247: DMGlobalToLocalEnd(stokes_da,X,INSERT_VALUES,X_local);
1248: DMDAVecGetArray(stokes_da,X_local,&stokes);
1250: DMDAGetInfo(stokes_da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
1251: DMDAGetBoundingBox(stokes_da,xymin,xymax);
1253: h = (xymax[0]-xymin[0])/((PetscReal)(M-1));
1255: tp_L2 = tu_L2 = tu_H1 = 0.0;
1256: tint_p_ms = tint_p = 0.0;
1258: DMDAGetElementCorners(stokes_da,&sex,&sey,&sez,&mx,&my,&mz);
1260: for (ek = sez; ek < sez+mz; ek++) {
1261: for (ej = sey; ej < sey+my; ej++) {
1262: for (ei = sex; ei < sex+mx; ei++) {
1263: /* get coords for the element */
1264: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1265: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1266: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1268: /* evaluate integral */
1269: for (p = 0; p < ngp; p++) {
1270: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1271: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1272: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1273: fac = gp_weight[p]*J_p;
1275: for (i = 0; i < NODES_PER_EL; i++) {
1276: tint_p_ms = tint_p_ms+fac*Ni_p[i]*stokes_analytic_e[i].p_dof;
1277: tint_p = tint_p +fac*Ni_p[i]*stokes_e[i].p_dof;
1278: }
1279: }
1280: }
1281: }
1282: }
1283: MPI_Allreduce(&tint_p_ms,&int_p_ms,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1284: MPI_Allreduce(&tint_p,&int_p,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1286: PetscPrintf(PETSC_COMM_WORLD,"\\int P dv %1.4e (h) %1.4e (ms)\n",PetscRealPart(int_p),PetscRealPart(int_p_ms));
1288: /* remove mine and add manufacture one */
1289: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1290: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1292: {
1293: PetscInt k,L,dof;
1294: PetscScalar *fields;
1295: DMDAGetInfo(stokes_da,0, 0,0,0, 0,0,0, &dof,0,0,0,0,0);
1297: VecGetLocalSize(X_local,&L);
1298: VecGetArray(X_local,&fields);
1299: for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1300: VecRestoreArray(X_local,&fields);
1302: VecGetLocalSize(X,&L);
1303: VecGetArray(X,&fields);
1304: for (k=0; k<L/dof; k++) fields[dof*k+3] = fields[dof*k+3] - int_p + int_p_ms;
1305: VecRestoreArray(X,&fields);
1306: }
1308: DMDAVecGetArray(stokes_da,X_local,&stokes);
1309: DMDAVecGetArray(stokes_da,X_analytic_local,&stokes_analytic);
1311: for (ek = sez; ek < sez+mz; ek++) {
1312: for (ej = sey; ej < sey+my; ej++) {
1313: for (ei = sex; ei < sex+mx; ei++) {
1314: /* get coords for the element */
1315: GetElementCoords3D(_coords,ei,ej,ek,el_coords);
1316: StokesDAGetNodalFields3D(stokes,ei,ej,ek,stokes_e);
1317: StokesDAGetNodalFields3D(stokes_analytic,ei,ej,ek,stokes_analytic_e);
1319: /* evaluate integral */
1320: p_e_L2 = 0.0;
1321: u_e_L2 = 0.0;
1322: u_e_H1 = 0.0;
1323: for (p = 0; p < ngp; p++) {
1324: ShapeFunctionQ13D_Evaluate(gp_xi[p],Ni_p);
1325: ShapeFunctionQ13D_Evaluate_dxi(gp_xi[p],GNi_p);
1326: ShapeFunctionQ13D_Evaluate_dx(GNi_p,GNx_p,el_coords,&J_p);
1327: fac = gp_weight[p]*J_p;
1329: for (i = 0; i < NODES_PER_EL; i++) {
1330: PetscScalar u_error,v_error,w_error;
1332: 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);
1334: u_error = stokes_e[i].u_dof-stokes_analytic_e[i].u_dof;
1335: v_error = stokes_e[i].v_dof-stokes_analytic_e[i].v_dof;
1336: w_error = stokes_e[i].w_dof-stokes_analytic_e[i].w_dof;
1337: /*
1338: if (p==0) {
1339: printf("p=0: %d %d %d %1.4e,%1.4e,%1.4e \n", ei,ej,ek,u_error,v_error,w_error);
1340: }
1341: */
1342: u_e_L2 += fac*Ni_p[i]*(u_error*u_error+v_error*v_error+w_error*w_error);
1344: u_e_H1 = u_e_H1+fac*(GNx_p[0][i]*u_error*GNx_p[0][i]*u_error /* du/dx */
1345: +GNx_p[1][i]*u_error*GNx_p[1][i]*u_error
1346: +GNx_p[2][i]*u_error*GNx_p[2][i]*u_error
1347: +GNx_p[0][i]*v_error*GNx_p[0][i]*v_error /* dv/dx */
1348: +GNx_p[1][i]*v_error*GNx_p[1][i]*v_error
1349: +GNx_p[2][i]*v_error*GNx_p[2][i]*v_error
1350: +GNx_p[0][i]*w_error*GNx_p[0][i]*w_error /* dw/dx */
1351: +GNx_p[1][i]*w_error*GNx_p[1][i]*w_error
1352: +GNx_p[2][i]*w_error*GNx_p[2][i]*w_error);
1353: }
1354: }
1356: tp_L2 += p_e_L2;
1357: tu_L2 += u_e_L2;
1358: tu_H1 += u_e_H1;
1359: }
1360: }
1361: }
1362: MPI_Allreduce(&tp_L2,&p_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1363: MPI_Allreduce(&tu_L2,&u_L2,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1364: MPI_Allreduce(&tu_H1,&u_H1,1,MPIU_SCALAR,MPIU_SUM,PETSC_COMM_WORLD);
1365: p_L2 = PetscSqrtScalar(p_L2);
1366: u_L2 = PetscSqrtScalar(u_L2);
1367: u_H1 = PetscSqrtScalar(u_H1);
1369: PetscPrintf(PETSC_COMM_WORLD,"%1.4e %1.4e %1.4e %1.4e \n",PetscRealPart(h),PetscRealPart(p_L2),PetscRealPart(u_L2),PetscRealPart(u_H1));
1372: DMDAVecRestoreArray(cda,coords,&_coords);
1374: DMDAVecRestoreArray(stokes_da,X_analytic_local,&stokes_analytic);
1375: VecDestroy(&X_analytic_local);
1376: DMDAVecRestoreArray(stokes_da,X_local,&stokes);
1377: VecDestroy(&X_local);
1378: return(0);
1379: }
1383: PetscErrorCode DAView_3DVTK_StructuredGrid_appended(DM da,Vec FIELD,const char file_prefix[])
1384: {
1385: char vtk_filename[PETSC_MAX_PATH_LEN];
1386: PetscMPIInt rank;
1387: MPI_Comm comm;
1388: FILE *vtk_fp = NULL;
1389: PetscInt si,sj,sk,nx,ny,nz,i;
1390: PetscInt f,n_fields,N;
1391: DM cda;
1392: Vec coords;
1393: Vec l_FIELD;
1394: PetscScalar *_L_FIELD;
1395: PetscInt memory_offset;
1396: PetscScalar *buffer;
1397: PetscLogDouble t0,t1;
1401: PetscTime(&t0);
1403: /* create file name */
1404: PetscObjectGetComm((PetscObject)da,&comm);
1405: MPI_Comm_rank(comm,&rank);
1406: PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"subdomain-%s-p%1.4d.vts",file_prefix,rank);
1408: /* open file and write header */
1409: vtk_fp = fopen(vtk_filename,"w");
1410: if (!vtk_fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename);
1412: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");
1414: /* coords */
1415: DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1416: N = nx * ny * nz;
1418: #if defined(PETSC_WORDS_BIGENDIAN)
1419: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1420: #else
1421: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"StructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1422: #endif
1423: 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);
1424: 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);
1426: memory_offset = 0;
1428: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <CellData></CellData>\n");
1430: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <Points>\n");
1432: /* copy coordinates */
1433: DMGetCoordinateDM(da,&cda);
1434: DMGetCoordinatesLocal(da,&coords);
1435: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"appended\" offset=\"%d\" />\n",memory_offset);
1436: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N*3;
1438: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </Points>\n");
1440: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PointData Scalars=\" ");
1441: DMDAGetInfo(da,0,0,0,0,0,0,0,&n_fields,0,0,0,0,0);
1442: for (f=0; f<n_fields; f++) {
1443: const char *field_name;
1444: DMDAGetFieldName(da,f,&field_name);
1445: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"%s ",field_name);
1446: }
1447: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\">\n");
1449: for (f=0; f<n_fields; f++) {
1450: const char *field_name;
1452: DMDAGetFieldName(da,f,&field_name);
1453: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <DataArray type=\"Float64\" Name=\"%s\" format=\"appended\" offset=\"%d\"/>\n", field_name,memory_offset);
1454: memory_offset = memory_offset + sizeof(PetscInt) + sizeof(PetscScalar)*N;
1455: }
1457: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PointData>\n");
1458: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </Piece>\n");
1459: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </StructuredGrid>\n");
1461: PetscMalloc(sizeof(PetscScalar)*N,&buffer);
1462: DMGetLocalVector(da,&l_FIELD);
1463: DMGlobalToLocalBegin(da, FIELD,INSERT_VALUES,l_FIELD);
1464: DMGlobalToLocalEnd(da,FIELD,INSERT_VALUES,l_FIELD);
1465: VecGetArray(l_FIELD,&_L_FIELD);
1467: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <AppendedData encoding=\"raw\">\n");
1468: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"_");
1470: /* write coordinates */
1471: {
1472: int length = sizeof(PetscScalar)*N*3;
1473: PetscScalar *allcoords;
1475: fwrite(&length,sizeof(int),1,vtk_fp);
1476: VecGetArray(coords,&allcoords);
1477: fwrite(allcoords,sizeof(PetscScalar),3*N,vtk_fp);
1478: VecRestoreArray(coords,&allcoords);
1479: }
1480: /* write fields */
1481: for (f=0; f<n_fields; f++) {
1482: int length = sizeof(PetscScalar)*N;
1483: fwrite(&length,sizeof(int),1,vtk_fp);
1484: /* load */
1485: for (i=0; i<N; i++) buffer[i] = _L_FIELD[n_fields*i + f];
1487: /* write */
1488: fwrite(buffer,sizeof(PetscScalar),N,vtk_fp);
1489: }
1490: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"\n </AppendedData>\n");
1492: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");
1494: PetscFree(buffer);
1495: VecRestoreArray(l_FIELD,&_L_FIELD);
1496: DMRestoreLocalVector(da,&l_FIELD);
1498: if (vtk_fp) {
1499: fclose(vtk_fp);
1500: vtk_fp = NULL;
1501: }
1503: PetscTime(&t1);
1504: return(0);
1505: }
1509: PetscErrorCode DAViewVTK_write_PieceExtend(FILE *vtk_fp,PetscInt indent_level,DM da,const char local_file_prefix[])
1510: {
1511: PetscMPIInt nproc,rank;
1512: MPI_Comm comm;
1513: const PetscInt *lx,*ly,*lz;
1514: PetscInt M,N,P,pM,pN,pP,sum,*olx,*oly,*olz;
1515: PetscInt *osx,*osy,*osz,*oex,*oey,*oez;
1516: PetscInt i,j,k,II,stencil;
1520: /* create file name */
1521: PetscObjectGetComm((PetscObject)da,&comm);
1522: MPI_Comm_size(comm,&nproc);
1523: MPI_Comm_rank(comm,&rank);
1525: DMDAGetInfo(da,0,&M,&N,&P,&pM,&pN,&pP,0,&stencil,0,0,0,0);
1526: DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
1528: /* generate start,end list */
1529: PetscMalloc(sizeof(PetscInt)*(pM+1),&olx);
1530: PetscMalloc(sizeof(PetscInt)*(pN+1),&oly);
1531: PetscMalloc(sizeof(PetscInt)*(pP+1),&olz);
1532: sum = 0;
1533: for (i=0; i<pM; i++) {
1534: olx[i] = sum;
1535: sum = sum + lx[i];
1536: }
1537: olx[pM] = sum;
1538: sum = 0;
1539: for (i=0; i<pN; i++) {
1540: oly[i] = sum;
1541: sum = sum + ly[i];
1542: }
1543: oly[pN] = sum;
1544: sum = 0;
1545: for (i=0; i<pP; i++) {
1546: olz[i] = sum;
1547: sum = sum + lz[i];
1548: }
1549: olz[pP] = sum;
1551: PetscMalloc(sizeof(PetscInt)*(pM),&osx);
1552: PetscMalloc(sizeof(PetscInt)*(pN),&osy);
1553: PetscMalloc(sizeof(PetscInt)*(pP),&osz);
1554: PetscMalloc(sizeof(PetscInt)*(pM),&oex);
1555: PetscMalloc(sizeof(PetscInt)*(pN),&oey);
1556: PetscMalloc(sizeof(PetscInt)*(pP),&oez);
1557: for (i=0; i<pM; i++) {
1558: osx[i] = olx[i] - stencil;
1559: oex[i] = olx[i] + lx[i] + stencil;
1560: if (osx[i]<0) osx[i]=0;
1561: if (oex[i]>M) oex[i]=M;
1562: }
1564: for (i=0; i<pN; i++) {
1565: osy[i] = oly[i] - stencil;
1566: oey[i] = oly[i] + ly[i] + stencil;
1567: if (osy[i]<0)osy[i]=0;
1568: if (oey[i]>M)oey[i]=N;
1569: }
1570: for (i=0; i<pP; i++) {
1571: osz[i] = olz[i] - stencil;
1572: oez[i] = olz[i] + lz[i] + stencil;
1573: if (osz[i]<0) osz[i]=0;
1574: if (oez[i]>P) oez[i]=P;
1575: }
1577: for (k=0; k<pP; k++) {
1578: for (j=0; j<pN; j++) {
1579: for (i=0; i<pM; i++) {
1580: char name[PETSC_MAX_PATH_LEN];
1581: PetscInt procid = i + j*pM + k*pM*pN; /* convert proc(i,j,k) to pid */
1582: PetscSNPrintf(name,sizeof(name),"subdomain-%s-p%1.4d.vts",local_file_prefix,procid);
1583: for (II=0; II<indent_level; II++) PetscFPrintf(PETSC_COMM_SELF,vtk_fp," ");
1585: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<Piece Extent=\"%d %d %d %d %d %d\" Source=\"%s\"/>\n",
1586: osx[i],oex[i]-1,
1587: osy[j],oey[j]-1,
1588: osz[k],oez[k]-1,name);
1589: }
1590: }
1591: }
1592: PetscFree(olx);
1593: PetscFree(oly);
1594: PetscFree(olz);
1595: PetscFree(osx);
1596: PetscFree(osy);
1597: PetscFree(osz);
1598: PetscFree(oex);
1599: PetscFree(oey);
1600: PetscFree(oez);
1601: return(0);
1602: }
1606: PetscErrorCode DAView_3DVTK_PStructuredGrid(DM da,const char file_prefix[],const char local_file_prefix[])
1607: {
1608: MPI_Comm comm;
1609: PetscMPIInt nproc,rank;
1610: char vtk_filename[PETSC_MAX_PATH_LEN];
1611: FILE *vtk_fp = NULL;
1612: PetscInt M,N,P,si,sj,sk,nx,ny,nz;
1613: PetscInt i,dofs;
1617: /* only master generates this file */
1618: PetscObjectGetComm((PetscObject)da,&comm);
1619: MPI_Comm_size(comm,&nproc);
1620: MPI_Comm_rank(comm,&rank);
1622: if (rank != 0) return(0);
1624: /* create file name */
1625: PetscSNPrintf(vtk_filename,sizeof(vtk_filename),"%s.pvts",file_prefix);
1626: vtk_fp = fopen(vtk_filename,"w");
1627: if (!vtk_fp) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SYS,"Cannot open file = %s \n",vtk_filename);
1629: /* (VTK) generate pvts header */
1630: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<?xml version=\"1.0\"?>\n");
1632: #if defined(PETSC_WORDS_BIGENDIAN)
1633: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"BigEndian\">\n");
1634: #else
1635: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"<VTKFile type=\"PStructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
1636: #endif
1638: /* define size of the nodal mesh based on the cell DM */
1639: DMDAGetInfo(da,0,&M,&N,&P,0,0,0,&dofs,0,0,0,0,0);
1640: DMDAGetGhostCorners(da,&si,&sj,&sk,&nx,&ny,&nz);
1641: 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 */
1643: /* DUMP THE CELL REFERENCES */
1644: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PCellData>\n");
1645: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PCellData>\n");
1647: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PPoints>\n");
1648: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"%d\"/>\n",NSD);
1649: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PPoints>\n");
1651: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PPointData>\n");
1652: for (i=0; i<dofs; i++) {
1653: const char *fieldname;
1654: DMDAGetFieldName(da,i,&fieldname);
1655: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," <PDataArray type=\"Float64\" Name=\"%s\" NumberOfComponents=\"1\"/>\n",fieldname);
1656: }
1657: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PPointData>\n");
1659: /* write out the parallel information */
1660: DAViewVTK_write_PieceExtend(vtk_fp,2,da,local_file_prefix);
1662: /* close the file */
1663: PetscFPrintf(PETSC_COMM_SELF,vtk_fp," </PStructuredGrid>\n");
1664: PetscFPrintf(PETSC_COMM_SELF,vtk_fp,"</VTKFile>\n");
1666: if (vtk_fp) {
1667: fclose(vtk_fp);
1668: vtk_fp = NULL;
1669: }
1670: return(0);
1671: }
1675: PetscErrorCode DAView3DPVTS(DM da, Vec x,const char NAME[])
1676: {
1677: char vts_filename[PETSC_MAX_PATH_LEN];
1678: char pvts_filename[PETSC_MAX_PATH_LEN];
1682: PetscSNPrintf(vts_filename,sizeof(vts_filename),"%s-mesh",NAME);
1683: DAView_3DVTK_StructuredGrid_appended(da,x,vts_filename);
1685: PetscSNPrintf(pvts_filename,sizeof(pvts_filename),"%s-mesh",NAME);
1686: DAView_3DVTK_PStructuredGrid(da,pvts_filename,vts_filename);
1687: return(0);
1688: }
1692: PetscErrorCode KSPMonitorStokesBlocks(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
1693: {
1695: PetscReal norms[4];
1696: Vec Br,v,w;
1697: Mat A;
1700: KSPGetOperators(ksp,&A,0,0);
1701: MatGetVecs(A,&w,&v);
1703: KSPBuildResidual(ksp,v,w,&Br);
1705: VecStrideNorm(Br,0,NORM_2,&norms[0]);
1706: VecStrideNorm(Br,1,NORM_2,&norms[1]);
1707: VecStrideNorm(Br,2,NORM_2,&norms[2]);
1708: VecStrideNorm(Br,3,NORM_2,&norms[3]);
1710: VecDestroy(&v);
1711: VecDestroy(&w);
1713: 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]);
1714: return(0);
1715: }
1719: static PetscErrorCode PCMGSetupViaCoarsen(PC pc,DM da_fine)
1720: {
1721: PetscInt nlevels,k,PETSC_UNUSED finest;
1722: DM *da_list,*daclist;
1723: Mat R;
1727: nlevels = 1;
1728: PetscOptionsGetInt(NULL,"-levels",&nlevels,0);
1730: PetscMalloc(sizeof(DM)*nlevels,&da_list);
1731: for (k=0; k<nlevels; k++) da_list[k] = NULL;
1732: PetscMalloc(sizeof(DM)*nlevels,&daclist);
1733: for (k=0; k<nlevels; k++) daclist[k] = NULL;
1735: /* finest grid is nlevels - 1 */
1736: finest = nlevels - 1;
1737: daclist[0] = da_fine;
1738: PetscObjectReference((PetscObject)da_fine);
1739: DMCoarsenHierarchy(da_fine,nlevels-1,&daclist[1]);
1740: for (k=0; k<nlevels; k++) {
1741: da_list[k] = daclist[nlevels-1-k];
1742: DMDASetUniformCoordinates(da_list[k],0.0,1.0,0.0,1.0,0.0,1.0);
1743: }
1745: PCMGSetLevels(pc,nlevels,NULL);
1746: PCMGSetType(pc,PC_MG_MULTIPLICATIVE);
1747: PCMGSetGalerkin(pc,PETSC_TRUE);
1749: for (k=1; k<nlevels; k++) {
1750: DMCreateInterpolation(da_list[k-1],da_list[k],&R,NULL);
1751: PCMGSetInterpolation(pc,k,R);
1752: MatDestroy(&R);
1753: }
1755: /* tidy up */
1756: for (k=0; k<nlevels; k++) {
1757: DMDestroy(&da_list[k]);
1758: }
1759: PetscFree(da_list);
1760: PetscFree(daclist);
1761: return(0);
1762: }
1766: static PetscErrorCode solve_stokes_3d_coupled(PetscInt mx,PetscInt my,PetscInt mz)
1767: {
1768: DM da_Stokes;
1769: PetscInt u_dof,p_dof,dof,stencil_width;
1770: Mat A,B;
1771: PetscInt mxl,myl,mzl;
1772: DM vel_cda;
1773: Vec vel_coords;
1774: PetscInt p;
1775: Vec f,X;
1776: DMDACoor3d ***_vel_coords;
1777: PetscInt its;
1778: KSP ksp_S;
1779: PetscInt model_definition = 0;
1780: PetscInt ei,ej,ek,sex,sey,sez,Imx,Jmy,Kmz;
1781: CellProperties cell_properties;
1782: PetscBool write_output = PETSC_FALSE;
1786: /* Generate the da for velocity and pressure */
1787: /* Num nodes in each direction is mx+1, my+1, mz+1 */
1788: u_dof = U_DOFS; /* Vx, Vy - velocities */
1789: p_dof = P_DOFS; /* p - pressure */
1790: dof = u_dof+p_dof;
1791: stencil_width = 1;
1792: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,
1793: mx+1,my+1,mz+1,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,NULL,&da_Stokes);
1794: DMDASetFieldName(da_Stokes,0,"Vx");
1795: DMDASetFieldName(da_Stokes,1,"Vy");
1796: DMDASetFieldName(da_Stokes,2,"Vz");
1797: DMDASetFieldName(da_Stokes,3,"P");
1799: /* unit box [0,1] x [0,1] x [0,1] */
1800: DMDASetUniformCoordinates(da_Stokes,0.0,1.0,0.0,1.0,0.0,1.0);
1802: /* local number of elements */
1803: DMDAGetLocalElementSize(da_Stokes,&mxl,&myl,&mzl);
1805: /* create quadrature point info for PDE coefficients */
1806: CellPropertiesCreate(da_Stokes,&cell_properties);
1808: /* interpolate the coordinates to quadrature points */
1809: DMGetCoordinateDM(da_Stokes,&vel_cda);
1810: DMGetCoordinatesLocal(da_Stokes,&vel_coords);
1811: DMDAVecGetArray(vel_cda,vel_coords,&_vel_coords);
1812: DMDAGetElementCorners(da_Stokes,&sex,&sey,&sez,&Imx,&Jmy,&Kmz);
1813: for (ek = sez; ek < sez+Kmz; ek++) {
1814: for (ej = sey; ej < sey+Jmy; ej++) {
1815: for (ei = sex; ei < sex+Imx; ei++) {
1816: /* get coords for the element */
1817: PetscInt ngp;
1818: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
1819: PetscScalar el_coords[NSD*NODES_PER_EL];
1820: GaussPointCoefficients *cell;
1822: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1823: GetElementCoords3D(_vel_coords,ei,ej,ek,el_coords);
1824: ConstructGaussQuadrature3D(&ngp,gp_xi,gp_weight);
1826: for (p = 0; p < GAUSS_POINTS; p++) {
1827: PetscScalar xi_p[NSD],Ni_p[NODES_PER_EL];
1828: PetscScalar gp_x,gp_y,gp_z;
1829: PetscInt n;
1831: xi_p[0] = gp_xi[p][0];
1832: xi_p[1] = gp_xi[p][1];
1833: xi_p[2] = gp_xi[p][2];
1834: ShapeFunctionQ13D_Evaluate(xi_p,Ni_p);
1836: gp_x = gp_y = gp_z = 0.0;
1837: for (n = 0; n < NODES_PER_EL; n++) {
1838: gp_x = gp_x+Ni_p[n]*el_coords[NSD*n];
1839: gp_y = gp_y+Ni_p[n]*el_coords[NSD*n+1];
1840: gp_z = gp_z+Ni_p[n]*el_coords[NSD*n+2];
1841: }
1842: cell->gp_coords[NSD*p] = gp_x;
1843: cell->gp_coords[NSD*p+1] = gp_y;
1844: cell->gp_coords[NSD*p+2] = gp_z;
1845: }
1846: }
1847: }
1848: }
1850: PetscOptionsGetInt(NULL,"-model",&model_definition,NULL);
1852: switch (model_definition) {
1853: case 0: /* isoviscous */
1854: for (ek = sez; ek < sez+Kmz; ek++) {
1855: for (ej = sey; ej < sey+Jmy; ej++) {
1856: for (ei = sex; ei < sex+Imx; ei++) {
1857: GaussPointCoefficients *cell;
1859: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1860: for (p = 0; p < GAUSS_POINTS; p++) {
1861: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1862: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1863: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1865: cell->eta[p] = 1.0;
1867: cell->fx[p] = 0.0*coord_x;
1868: cell->fy[p] = -sin((double)2.2*M_PI*coord_y)*cos(1.0*M_PI*coord_x);
1869: cell->fz[p] = 0.0*coord_z;
1870: cell->hc[p] = 0.0;
1871: }
1872: }
1873: }
1874: }
1875: break;
1877: case 1: /* manufactured */
1878: for (ek = sez; ek < sez+Kmz; ek++) {
1879: for (ej = sey; ej < sey+Jmy; ej++) {
1880: for (ei = sex; ei < sex+Imx; ei++) {
1881: PetscReal eta,Fm[NSD],Fc,pos[NSD];
1882: GaussPointCoefficients *cell;
1884: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1885: for (p = 0; p < GAUSS_POINTS; p++) {
1886: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1887: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1888: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1890: pos[0] = coord_x;
1891: pos[1] = coord_y;
1892: pos[2] = coord_z;
1894: evaluate_MS_FrankKamentski(pos,NULL,NULL,&eta,Fm,&Fc);
1895: cell->eta[p] = eta;
1897: cell->fx[p] = Fm[0];
1898: cell->fy[p] = Fm[1];
1899: cell->fz[p] = Fm[2];
1900: cell->hc[p] = Fc;
1901: }
1902: }
1903: }
1904: }
1905: break;
1907: case 2: /* solcx */
1908: for (ek = sez; ek < sez+Kmz; ek++) {
1909: for (ej = sey; ej < sey+Jmy; ej++) {
1910: for (ei = sex; ei < sex+Imx; ei++) {
1911: GaussPointCoefficients *cell;
1913: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1914: for (p = 0; p < GAUSS_POINTS; p++) {
1915: PetscReal coord_x = PetscRealPart(cell->gp_coords[NSD*p]);
1916: PetscReal coord_y = PetscRealPart(cell->gp_coords[NSD*p+1]);
1917: PetscReal coord_z = PetscRealPart(cell->gp_coords[NSD*p+2]);
1919: cell->eta[p] = 1.0;
1921: cell->fx[p] = 0.0;
1922: cell->fy[p] = -sin((double)3*M_PI*coord_y)*cos(1.0*M_PI*coord_x);
1923: cell->fz[p] = 0.0*coord_z;
1924: cell->hc[p] = 0.0;
1925: }
1926: }
1927: }
1928: }
1929: break;
1931: case 3: /* sinker */
1932: for (ek = sez; ek < sez+Kmz; ek++) {
1933: for (ej = sey; ej < sey+Jmy; ej++) {
1934: for (ei = sex; ei < sex+Imx; ei++) {
1935: GaussPointCoefficients *cell;
1937: CellPropertiesGetCell(cell_properties,ei,ej,ek,&cell);
1938: for (p = 0; p < GAUSS_POINTS; p++) {
1939: PetscReal xp = PetscRealPart(cell->gp_coords[NSD*p]);
1940: PetscReal yp = PetscRealPart(cell->gp_coords[NSD*p+1]);
1941: PetscReal zp = PetscRealPart(cell->gp_coords[NSD*p+2]);
1943: cell->eta[p] = 1.0e-2;
1944: cell->fx[p] = 0.0;
1945: cell->fy[p] = 0.0;
1946: cell->fz[p] = 0.0;
1947: cell->hc[p] = 0.0;
1949: if ((fabs(xp-0.5) < 0.2) &&
1950: (fabs(yp-0.5) < 0.2) &&
1951: (fabs(zp-0.5) < 0.2)) {
1952: cell->eta[p] = 1.0;
1953: cell->fz[p] = 1.0;
1954: }
1956: }
1957: }
1958: }
1959: }
1960: break;
1962: default:
1963: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"No default model is supported. Choose either -model {0,1,2,3}");
1964: break;
1965: }
1967: DMDAVecRestoreArray(vel_cda,vel_coords,&_vel_coords);
1969: /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1970: DMCreateMatrix(da_Stokes,MATAIJ,&A);
1971: DMCreateMatrix(da_Stokes,MATAIJ,&B);
1972: DMCreateGlobalVector(da_Stokes,&X);
1973: DMCreateGlobalVector(da_Stokes,&f);
1975: /* assemble A11 */
1976: MatZeroEntries(A);
1977: MatZeroEntries(B);
1978: VecZeroEntries(f);
1980: AssembleA_Stokes(A,da_Stokes,cell_properties);
1981: AssembleA_PCStokes(B,da_Stokes,cell_properties);
1982: /* build force vector */
1983: AssembleF_Stokes(f,da_Stokes,cell_properties);
1985: /* SOLVE */
1986: KSPCreate(PETSC_COMM_WORLD,&ksp_S);
1987: KSPSetOptionsPrefix(ksp_S,"stokes_"); /* stokes */
1988: KSPSetOperators(ksp_S,A,B,SAME_NONZERO_PATTERN);
1989: KSPSetFromOptions(ksp_S);
1991: {
1992: PC pc;
1993: const PetscInt ufields[] = {0,1,2},pfields[] = {3};
1994: KSPGetPC(ksp_S,&pc);
1995: PCFieldSplitSetBlockSize(pc,4);
1996: PCFieldSplitSetFields(pc,"u",3,ufields,ufields);
1997: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1998: }
2000: {
2001: PC pc;
2002: PetscBool same = PETSC_FALSE;
2003: KSPGetPC(ksp_S,&pc);
2004: PetscObjectTypeCompare((PetscObject)pc,PCMG,&same);
2005: if (same) {
2006: PCMGSetupViaCoarsen(pc,da_Stokes);
2007: }
2008: }
2010: {
2011: PetscBool stokes_monitor = PETSC_FALSE;
2012: PetscOptionsGetBool(NULL,"-stokes_ksp_monitor_blocks",&stokes_monitor,0);
2013: if (stokes_monitor) {
2014: KSPMonitorSet(ksp_S,KSPMonitorStokesBlocks,NULL,NULL);
2015: }
2016: }
2017: KSPSolve(ksp_S,f,X);
2019: PetscOptionsGetBool(NULL,"-write_pvts",&write_output,NULL);
2020: if (write_output) {DAView3DPVTS(da_Stokes,X,"up");}
2021: {
2022: PetscBool flg = PETSC_FALSE;
2023: char filename[PETSC_MAX_PATH_LEN];
2024: PetscOptionsGetString(NULL,"-write_binary",filename,sizeof(filename),&flg);
2025: if (flg) {
2026: PetscViewer viewer;
2027: /* PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename[0]?filename:"ex42-binaryoutput",FILE_MODE_WRITE,&viewer); */
2028: PetscViewerVTKOpen(PETSC_COMM_WORLD,"ex42.vts",FILE_MODE_WRITE,&viewer);
2029: VecView(X,viewer);
2030: PetscViewerDestroy(&viewer);
2031: }
2032: }
2033: KSPGetIterationNumber(ksp_S,&its);
2035: /* verify */
2036: if (model_definition == 1) {
2037: DM da_Stokes_analytic;
2038: Vec X_analytic;
2040: DMDACreateManufacturedSolution(mx,my,mz,&da_Stokes_analytic,&X_analytic);
2041: if (write_output) {
2042: DAView3DPVTS(da_Stokes_analytic,X_analytic,"ms");
2043: }
2044: DMDAIntegrateErrors3D(da_Stokes_analytic,X,X_analytic);
2045: if (write_output) {
2046: DAView3DPVTS(da_Stokes,X,"up2");
2047: }
2048: DMDestroy(&da_Stokes_analytic);
2049: VecDestroy(&X_analytic);
2050: }
2052: KSPDestroy(&ksp_S);
2053: VecDestroy(&X);
2054: VecDestroy(&f);
2055: MatDestroy(&A);
2056: MatDestroy(&B);
2058: CellPropertiesDestroy(&cell_properties);
2059: DMDestroy(&da_Stokes);
2060: return(0);
2061: }
2065: int main(int argc,char **args)
2066: {
2068: PetscInt mx,my,mz;
2070: PetscInitialize(&argc,&args,(char*)0,help);
2072: mx = my = mz = 10;
2073: PetscOptionsGetInt(NULL,"-mx",&mx,NULL);
2074: my = mx; mz = mx;
2075: PetscOptionsGetInt(NULL,"-my",&my,NULL);
2076: PetscOptionsGetInt(NULL,"-mz",&mz,NULL);
2078: solve_stokes_3d_coupled(mx,my,mz);
2080: PetscFinalize();
2081: return 0;
2082: }