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