1 static char help[] = "Solves the incompressible, variable viscosity Stokes equation in 2d on the unit domain \n\
2 using Q1Q1 elements, stabilized with Bochev's polynomial projection method. \n\
3 The models defined utilise free slip boundary conditions on all sides. \n\
4 Options: \n"
5 "\
6 -mx : Number of elements in the x-direction \n\
7 -my : Number of elements in the y-direction \n\
8 -o : Specify output filename for solution (will be PETSc binary format or paraview format if the extension is .vts) \n\
9 -gnuplot : Output Gauss point coordinates, coefficients and u,p solution in gnuplot format \n\
10 -glvis : Visualizes coefficients and u,p solution through GLVIs (use -viewer_glvis_dmda_bs 2,1 to visualize velocity as a vector)\n\
11 -c_str : Indicates the structure of the coefficients to use \n"
12 "\
13 -c_str 0 => Coefficient definition for an analytic solution with a vertical jump in viscosity at x = xc \n\
14 This problem is driven by the forcing function f(x,y) = (0, sin(nz pi y)cos(pi x) \n\
15 Parameters: \n\
16 -solcx_eta0 : Viscosity to the left of the interface \n\
17 -solcx_eta1 : Viscosity to the right of the interface \n\
18 -solcx_xc : Location of the interface \n\
19 -solcx_nz : Wavenumber in the y direction \n"
20 "\
21 -c_str 1 => Coefficient definition for a dense rectangular blob located at the center of the domain \n\
22 Parameters: \n\
23 -sinker_eta0 : Viscosity of the background fluid \n\
24 -sinker_eta1 : Viscosity of the blob \n\
25 -sinker_dx : Width of the blob \n\
26 -sinker_dy : Height of the blob \n"
27 "\
28 -c_str 2 => Coefficient definition for a dense circular blob located at the center of the domain \n\
29 Parameters: \n\
30 -sinker_eta0 : Viscosity of the background fluid \n\
31 -sinker_eta1 : Viscosity of the blob \n\
32 -sinker_r : Radius of the blob \n"
33 "\
34 -c_str 3 => Coefficient definition for a dense circular and rectangular inclusion (located at the center of the domain) \n\
35 -sinker_eta0 : Viscosity of the background fluid \n\
36 -sinker_eta1 : Viscosity of the two inclusions \n\
37 -sinker_r : Radius of the circular inclusion \n\
38 -sinker_c0x : Origin (x-coord) of the circular inclusion \n\
39 -sinker_c0y : Origin (y-coord) of the circular inclusion \n\
40 -sinker_dx : Width of the rectangular inclusion \n\
41 -sinker_dy : Height of the rectangular inclusion \n\
42 -sinker_phi : Rotation angle of the rectangular inclusion \n"
43 "\
44 -c_str 4 => Coefficient definition for checkerboard jumps aligned with the domain decomposition \n\
45 -jump_eta0 : Viscosity for black subdomains \n\
46 -jump_magnitude : Magnitude of jumps. White subdomains will have eta = eta0*10^magnitude \n\
47 -jump_nz : Wavenumber in the y direction for rhs \n"
48 "\
49 -use_gp_coords : Evaluate the viscosity and force term at the global coordinates of each quadrature point \n\
50 By default, the viscosity and force term are evaluated at the element center and applied as a constant over the entire element \n";
51
52 /* Contributed by Dave May */
53
54 #include <petscksp.h>
55 #include <petscdm.h>
56 #include <petscdmda.h>
57
58 /* A Maple-generated exact solution created by Mirko Velic (mirko.velic@sci.monash.edu.au) */
59 #include "ex43-solcx.h"
60
61 static PetscErrorCode DMDABCApplyFreeSlip(DM, Mat, Vec);
62
63 #define NSD 2 /* number of spatial dimensions */
64 #define NODES_PER_EL 4 /* nodes per element */
65 #define U_DOFS 2 /* degrees of freedom per velocity node */
66 #define P_DOFS 1 /* degrees of freedom per pressure node */
67 #define GAUSS_POINTS 4
68
69 /* Gauss point based evaluation 8+4+4+4 = 20 */
70 typedef struct {
71 PetscScalar gp_coords[2 * GAUSS_POINTS];
72 PetscScalar eta[GAUSS_POINTS];
73 PetscScalar fx[GAUSS_POINTS];
74 PetscScalar fy[GAUSS_POINTS];
75 } GaussPointCoefficients;
76
77 typedef struct {
78 PetscScalar u_dof;
79 PetscScalar v_dof;
80 PetscScalar p_dof;
81 } StokesDOF;
82
glvis_extract_eta(PetscObject oV,PetscInt nf,PetscObject oVf[],PetscCtx ctx)83 static PetscErrorCode glvis_extract_eta(PetscObject oV, PetscInt nf, PetscObject oVf[], PetscCtx ctx)
84 {
85 DM properties_da = (DM)ctx, stokes_da;
86 Vec V = (Vec)oV, *Vf = (Vec *)oVf;
87 GaussPointCoefficients **props;
88 PetscInt sex, sey, mx, my;
89 PetscInt ei, ej, p, cum;
90 PetscScalar *array;
91
92 PetscFunctionBeginUser;
93 PetscCall(VecGetDM(Vf[0], &stokes_da));
94 PetscCall(DMDAVecGetArrayRead(properties_da, V, &props));
95 PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
96 PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
97 PetscCall(VecGetArray(Vf[0], &array));
98 cum = 0;
99 for (ej = sey; ej < sey + my; ej++) {
100 for (ei = sex; ei < sex + mx; ei++) {
101 for (p = 0; p < GAUSS_POINTS; p++) array[cum++] = props[ej][ei].eta[p];
102 }
103 }
104 PetscCall(VecRestoreArray(Vf[0], &array));
105 PetscCall(DMDAVecRestoreArrayRead(properties_da, V, &props));
106 PetscFunctionReturn(PETSC_SUCCESS);
107 }
108
109 /*
110 Element: Local basis function ordering
111 1-----2
112 | |
113 | |
114 0-----3
115 */
ConstructQ12D_Ni(PetscScalar _xi[],PetscScalar Ni[])116 static void ConstructQ12D_Ni(PetscScalar _xi[], PetscScalar Ni[])
117 {
118 PetscScalar xi = _xi[0];
119 PetscScalar eta = _xi[1];
120
121 Ni[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
122 Ni[1] = 0.25 * (1.0 - xi) * (1.0 + eta);
123 Ni[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
124 Ni[3] = 0.25 * (1.0 + xi) * (1.0 - eta);
125 }
126
ConstructQ12D_GNi(PetscScalar _xi[],PetscScalar GNi[][NODES_PER_EL])127 static void ConstructQ12D_GNi(PetscScalar _xi[], PetscScalar GNi[][NODES_PER_EL])
128 {
129 PetscScalar xi = _xi[0];
130 PetscScalar eta = _xi[1];
131
132 GNi[0][0] = -0.25 * (1.0 - eta);
133 GNi[0][1] = -0.25 * (1.0 + eta);
134 GNi[0][2] = 0.25 * (1.0 + eta);
135 GNi[0][3] = 0.25 * (1.0 - eta);
136
137 GNi[1][0] = -0.25 * (1.0 - xi);
138 GNi[1][1] = 0.25 * (1.0 - xi);
139 GNi[1][2] = 0.25 * (1.0 + xi);
140 GNi[1][3] = -0.25 * (1.0 + xi);
141 }
142
ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL],PetscScalar GNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar * det_J)143 static void ConstructQ12D_GNx(PetscScalar GNi[][NODES_PER_EL], PetscScalar GNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
144 {
145 PetscScalar J00, J01, J10, J11, J;
146 PetscScalar iJ00, iJ01, iJ10, iJ11;
147 PetscInt i;
148
149 J00 = J01 = J10 = J11 = 0.0;
150 for (i = 0; i < NODES_PER_EL; i++) {
151 PetscScalar cx = coords[2 * i];
152 PetscScalar cy = coords[2 * i + 1];
153
154 J00 = J00 + GNi[0][i] * cx; /* J_xx = dx/dxi */
155 J01 = J01 + GNi[0][i] * cy; /* J_xy = dy/dxi */
156 J10 = J10 + GNi[1][i] * cx; /* J_yx = dx/deta */
157 J11 = J11 + GNi[1][i] * cy; /* J_yy = dy/deta */
158 }
159 J = (J00 * J11) - (J01 * J10);
160
161 iJ00 = J11 / J;
162 iJ01 = -J01 / J;
163 iJ10 = -J10 / J;
164 iJ11 = J00 / J;
165
166 for (i = 0; i < NODES_PER_EL; i++) {
167 GNx[0][i] = GNi[0][i] * iJ00 + GNi[1][i] * iJ01;
168 GNx[1][i] = GNi[0][i] * iJ10 + GNi[1][i] * iJ11;
169 }
170
171 if (det_J) *det_J = J;
172 }
173
ConstructGaussQuadrature(PetscInt * ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])174 static void ConstructGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
175 {
176 *ngp = 4;
177 gp_xi[0][0] = -0.57735026919;
178 gp_xi[0][1] = -0.57735026919;
179 gp_xi[1][0] = -0.57735026919;
180 gp_xi[1][1] = 0.57735026919;
181 gp_xi[2][0] = 0.57735026919;
182 gp_xi[2][1] = 0.57735026919;
183 gp_xi[3][0] = 0.57735026919;
184 gp_xi[3][1] = -0.57735026919;
185 gp_weight[0] = 1.0;
186 gp_weight[1] = 1.0;
187 gp_weight[2] = 1.0;
188 gp_weight[3] = 1.0;
189 }
190
191 /*
192 i,j are the element indices
193 The unknown is a vector quantity.
194 The s[].c is used to indicate the degree of freedom.
195 */
DMDAGetElementEqnums_up(MatStencil s_u[],MatStencil s_p[],PetscInt i,PetscInt j)196 static PetscErrorCode DMDAGetElementEqnums_up(MatStencil s_u[], MatStencil s_p[], PetscInt i, PetscInt j)
197 {
198 PetscFunctionBeginUser;
199 /* velocity */
200 /* node 0 */
201 s_u[0].i = i;
202 s_u[0].j = j;
203 s_u[0].c = 0; /* Vx0 */
204 s_u[1].i = i;
205 s_u[1].j = j;
206 s_u[1].c = 1; /* Vy0 */
207
208 /* node 1 */
209 s_u[2].i = i;
210 s_u[2].j = j + 1;
211 s_u[2].c = 0; /* Vx1 */
212 s_u[3].i = i;
213 s_u[3].j = j + 1;
214 s_u[3].c = 1; /* Vy1 */
215
216 /* node 2 */
217 s_u[4].i = i + 1;
218 s_u[4].j = j + 1;
219 s_u[4].c = 0; /* Vx2 */
220 s_u[5].i = i + 1;
221 s_u[5].j = j + 1;
222 s_u[5].c = 1; /* Vy2 */
223
224 /* node 3 */
225 s_u[6].i = i + 1;
226 s_u[6].j = j;
227 s_u[6].c = 0; /* Vx3 */
228 s_u[7].i = i + 1;
229 s_u[7].j = j;
230 s_u[7].c = 1; /* Vy3 */
231
232 /* pressure */
233 s_p[0].i = i;
234 s_p[0].j = j;
235 s_p[0].c = 2; /* P0 */
236 s_p[1].i = i;
237 s_p[1].j = j + 1;
238 s_p[1].c = 2; /* P1 */
239 s_p[2].i = i + 1;
240 s_p[2].j = j + 1;
241 s_p[2].c = 2; /* P2 */
242 s_p[3].i = i + 1;
243 s_p[3].j = j;
244 s_p[3].c = 2; /* P3 */
245 PetscFunctionReturn(PETSC_SUCCESS);
246 }
247
DMDAGetElementOwnershipRanges2d(DM da,PetscInt ** _lx,PetscInt ** _ly)248 static PetscErrorCode DMDAGetElementOwnershipRanges2d(DM da, PetscInt **_lx, PetscInt **_ly)
249 {
250 PetscMPIInt rank;
251 PetscInt proc_I, proc_J;
252 PetscInt cpu_x, cpu_y;
253 PetscInt local_mx, local_my;
254 Vec vlx, vly;
255 PetscInt *LX, *LY, i;
256 PetscScalar *_a;
257 Vec V_SEQ;
258 VecScatter ctx;
259
260 PetscFunctionBeginUser;
261 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
262
263 PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
264
265 proc_J = rank / cpu_x;
266 proc_I = rank - cpu_x * proc_J;
267
268 PetscCall(PetscMalloc1(cpu_x, &LX));
269 PetscCall(PetscMalloc1(cpu_y, &LY));
270
271 PetscCall(DMDAGetElementsSizes(da, &local_mx, &local_my, NULL));
272 PetscCall(VecCreate(PETSC_COMM_WORLD, &vlx));
273 PetscCall(VecSetSizes(vlx, PETSC_DECIDE, cpu_x));
274 PetscCall(VecSetFromOptions(vlx));
275
276 PetscCall(VecCreate(PETSC_COMM_WORLD, &vly));
277 PetscCall(VecSetSizes(vly, PETSC_DECIDE, cpu_y));
278 PetscCall(VecSetFromOptions(vly));
279
280 PetscCall(VecSetValue(vlx, proc_I, (PetscScalar)(local_mx + 1.0e-9), INSERT_VALUES));
281 PetscCall(VecSetValue(vly, proc_J, (PetscScalar)(local_my + 1.0e-9), INSERT_VALUES));
282 PetscCall(VecAssemblyBegin(vlx));
283 PetscCall(VecAssemblyEnd(vlx));
284 PetscCall(VecAssemblyBegin(vly));
285 PetscCall(VecAssemblyEnd(vly));
286
287 PetscCall(VecScatterCreateToAll(vlx, &ctx, &V_SEQ));
288 PetscCall(VecScatterBegin(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
289 PetscCall(VecScatterEnd(ctx, vlx, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
290 PetscCall(VecGetArray(V_SEQ, &_a));
291 for (i = 0; i < cpu_x; i++) LX[i] = (PetscInt)PetscRealPart(_a[i]);
292 PetscCall(VecRestoreArray(V_SEQ, &_a));
293 PetscCall(VecScatterDestroy(&ctx));
294 PetscCall(VecDestroy(&V_SEQ));
295
296 PetscCall(VecScatterCreateToAll(vly, &ctx, &V_SEQ));
297 PetscCall(VecScatterBegin(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
298 PetscCall(VecScatterEnd(ctx, vly, V_SEQ, INSERT_VALUES, SCATTER_FORWARD));
299 PetscCall(VecGetArray(V_SEQ, &_a));
300 for (i = 0; i < cpu_y; i++) LY[i] = (PetscInt)PetscRealPart(_a[i]);
301 PetscCall(VecRestoreArray(V_SEQ, &_a));
302 PetscCall(VecScatterDestroy(&ctx));
303 PetscCall(VecDestroy(&V_SEQ));
304
305 *_lx = LX;
306 *_ly = LY;
307
308 PetscCall(VecDestroy(&vlx));
309 PetscCall(VecDestroy(&vly));
310 PetscFunctionReturn(PETSC_SUCCESS);
311 }
312
DMDACoordViewGnuplot2d(DM da,const char prefix[])313 static PetscErrorCode DMDACoordViewGnuplot2d(DM da, const char prefix[])
314 {
315 DM cda;
316 Vec coords;
317 DMDACoor2d **_coords;
318 PetscInt si, sj, nx, ny, i, j;
319 FILE *fp;
320 char fname[PETSC_MAX_PATH_LEN];
321 PetscMPIInt rank;
322
323 PetscFunctionBeginUser;
324 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
325 PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
326 PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
327 PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
328
329 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### Element geometry for processor %1.4d ### \n", rank));
330
331 PetscCall(DMGetCoordinateDM(da, &cda));
332 PetscCall(DMGetCoordinatesLocal(da, &coords));
333 PetscCall(DMDAVecGetArrayRead(cda, coords, &_coords));
334 PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
335 for (j = sj; j < sj + ny - 1; j++) {
336 for (i = si; i < si + nx - 1; i++) {
337 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
338 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i].x), (double)PetscRealPart(_coords[j + 1][i].y)));
339 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j + 1][i + 1].x), (double)PetscRealPart(_coords[j + 1][i + 1].y)));
340 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n", (double)PetscRealPart(_coords[j][i + 1].x), (double)PetscRealPart(_coords[j][i + 1].y)));
341 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e \n\n", (double)PetscRealPart(_coords[j][i].x), (double)PetscRealPart(_coords[j][i].y)));
342 }
343 }
344 PetscCall(DMDAVecRestoreArrayRead(cda, coords, &_coords));
345
346 PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
347 PetscFunctionReturn(PETSC_SUCCESS);
348 }
349
DMDAViewGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])350 static PetscErrorCode DMDAViewGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
351 {
352 DM cda;
353 Vec coords, local_fields;
354 DMDACoor2d **_coords;
355 FILE *fp;
356 char fname[PETSC_MAX_PATH_LEN];
357 PetscMPIInt rank;
358 PetscInt si, sj, nx, ny, i, j;
359 PetscInt n_dofs, d;
360 PetscScalar *_fields;
361
362 PetscFunctionBeginUser;
363 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
364 PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
365 PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
366 PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
367
368 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
369 PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
370 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
371 for (d = 0; d < n_dofs; d++) {
372 const char *field_name;
373 PetscCall(DMDAGetFieldName(da, d, &field_name));
374 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
375 }
376 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));
377
378 PetscCall(DMGetCoordinateDM(da, &cda));
379 PetscCall(DMGetCoordinatesLocal(da, &coords));
380 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
381 PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
382
383 PetscCall(DMCreateLocalVector(da, &local_fields));
384 PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
385 PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
386 PetscCall(VecGetArray(local_fields, &_fields));
387
388 for (j = sj; j < sj + ny; j++) {
389 for (i = si; i < si + nx; i++) {
390 PetscScalar coord_x, coord_y;
391 PetscScalar field_d;
392
393 coord_x = _coords[j][i].x;
394 coord_y = _coords[j][i].y;
395
396 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
397 for (d = 0; d < n_dofs; d++) {
398 field_d = _fields[n_dofs * ((i - si) + (j - sj) * (nx)) + d];
399 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e ", (double)PetscRealPart(field_d)));
400 }
401 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
402 }
403 }
404 PetscCall(VecRestoreArray(local_fields, &_fields));
405 PetscCall(VecDestroy(&local_fields));
406
407 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
408
409 PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
410 PetscFunctionReturn(PETSC_SUCCESS);
411 }
412
DMDAViewCoefficientsGnuplot2d(DM da,Vec fields,const char comment[],const char prefix[])413 static PetscErrorCode DMDAViewCoefficientsGnuplot2d(DM da, Vec fields, const char comment[], const char prefix[])
414 {
415 DM cda;
416 Vec local_fields;
417 FILE *fp;
418 char fname[PETSC_MAX_PATH_LEN];
419 PetscMPIInt rank;
420 PetscInt si, sj, nx, ny, i, j, p;
421 PetscInt n_dofs, d;
422 GaussPointCoefficients **_coefficients;
423
424 PetscFunctionBeginUser;
425 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
426 PetscCall(PetscSNPrintf(fname, sizeof(fname), "%s-p%1.4d.dat", prefix, rank));
427 PetscCall(PetscFOpen(PETSC_COMM_SELF, fname, "w", &fp));
428 PetscCheck(fp, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot open file");
429
430 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### %s (processor %1.4d) ### \n", comment, rank));
431 PetscCall(DMDAGetInfo(da, 0, 0, 0, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
432 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "### x y "));
433 for (d = 0; d < n_dofs; d++) {
434 const char *field_name;
435 PetscCall(DMDAGetFieldName(da, d, &field_name));
436 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%s ", field_name));
437 }
438 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "###\n"));
439
440 PetscCall(DMGetCoordinateDM(da, &cda));
441 PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
442
443 PetscCall(DMCreateLocalVector(da, &local_fields));
444 PetscCall(DMGlobalToLocalBegin(da, fields, INSERT_VALUES, local_fields));
445 PetscCall(DMGlobalToLocalEnd(da, fields, INSERT_VALUES, local_fields));
446 PetscCall(DMDAVecGetArray(da, local_fields, &_coefficients));
447
448 for (j = sj; j < sj + ny; j++) {
449 for (i = si; i < si + nx; i++) {
450 PetscScalar coord_x, coord_y;
451
452 for (p = 0; p < GAUSS_POINTS; p++) {
453 coord_x = _coefficients[j][i].gp_coords[2 * p];
454 coord_y = _coefficients[j][i].gp_coords[2 * p + 1];
455
456 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e ", (double)PetscRealPart(coord_x), (double)PetscRealPart(coord_y)));
457
458 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "%1.6e %1.6e %1.6e", (double)PetscRealPart(_coefficients[j][i].eta[p]), (double)PetscRealPart(_coefficients[j][i].fx[p]), (double)PetscRealPart(_coefficients[j][i].fy[p])));
459 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fp, "\n"));
460 }
461 }
462 }
463 PetscCall(DMDAVecRestoreArray(da, local_fields, &_coefficients));
464 PetscCall(VecDestroy(&local_fields));
465
466 PetscCall(PetscFClose(PETSC_COMM_SELF, fp));
467 PetscFunctionReturn(PETSC_SUCCESS);
468 }
469
ASS_MAP_wIwDI_uJuDJ(PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)470 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)
471 {
472 PetscInt ij;
473 PetscInt r, c, nc;
474
475 nc = u_NPE * u_dof;
476 r = w_dof * wi + wd;
477 c = u_dof * ui + ud;
478 ij = r * nc + c;
479 return ij;
480 }
481
482 /*
483 D = [ 2.eta 0 0 ]
484 [ 0 2.eta 0 ]
485 [ 0 0 eta ]
486
487 B = [ d_dx 0 ]
488 [ 0 d_dy ]
489 [ d_dy d_dx ]
490 */
FormStressOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])491 static void FormStressOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
492 {
493 PetscInt ngp;
494 PetscScalar gp_xi[GAUSS_POINTS][2];
495 PetscScalar gp_weight[GAUSS_POINTS];
496 PetscInt p, i, j, k;
497 PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
498 PetscScalar J_p, tildeD[3];
499 PetscScalar B[3][U_DOFS * NODES_PER_EL];
500
501 /* define quadrature rule */
502 ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
503
504 /* evaluate integral */
505 for (p = 0; p < ngp; p++) {
506 ConstructQ12D_GNi(gp_xi[p], GNi_p);
507 ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
508
509 for (i = 0; i < NODES_PER_EL; i++) {
510 PetscScalar d_dx_i = GNx_p[0][i];
511 PetscScalar d_dy_i = GNx_p[1][i];
512
513 B[0][2 * i] = d_dx_i;
514 B[0][2 * i + 1] = 0.0;
515 B[1][2 * i] = 0.0;
516 B[1][2 * i + 1] = d_dy_i;
517 B[2][2 * i] = d_dy_i;
518 B[2][2 * i + 1] = d_dx_i;
519 }
520
521 tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
522 tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
523 tildeD[2] = gp_weight[p] * J_p * eta[p];
524
525 /* form Bt tildeD B */
526 /*
527 Ke_ij = Bt_ik . D_kl . B_lj
528 = B_ki . D_kl . B_lj
529 = B_ki . D_kk . B_kj
530 */
531 for (i = 0; i < 8; i++) {
532 for (j = 0; j < 8; j++) {
533 for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
534 Ke[i + 8 * j] = Ke[i + 8 * j] + B[k][i] * tildeD[k] * B[k][j];
535 }
536 }
537 }
538 }
539 }
540
FormGradientOperatorQ1(PetscScalar Ke[],PetscScalar coords[])541 static void FormGradientOperatorQ1(PetscScalar Ke[], PetscScalar coords[])
542 {
543 PetscInt ngp;
544 PetscScalar gp_xi[GAUSS_POINTS][2];
545 PetscScalar gp_weight[GAUSS_POINTS];
546 PetscInt p, i, j, di;
547 PetscScalar Ni_p[NODES_PER_EL];
548 PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
549 PetscScalar J_p, fac;
550
551 /* define quadrature rule */
552 ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
553
554 /* evaluate integral */
555 for (p = 0; p < ngp; p++) {
556 ConstructQ12D_Ni(gp_xi[p], Ni_p);
557 ConstructQ12D_GNi(gp_xi[p], GNi_p);
558 ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
559 fac = gp_weight[p] * J_p;
560
561 for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
562 for (di = 0; di < NSD; di++) { /* u dofs */
563 for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
564 PetscInt IJ;
565 IJ = ASS_MAP_wIwDI_uJuDJ(i, di, NODES_PER_EL, 2, j, 0, NODES_PER_EL, 1);
566
567 Ke[IJ] = Ke[IJ] - GNx_p[di][i] * Ni_p[j] * fac;
568 }
569 }
570 }
571 }
572 }
573
FormDivergenceOperatorQ1(PetscScalar De[],PetscScalar coords[])574 static void FormDivergenceOperatorQ1(PetscScalar De[], PetscScalar coords[])
575 {
576 PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
577 PetscInt i, j;
578 PetscInt nr_g, nc_g;
579
580 PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));
581 FormGradientOperatorQ1(Ge, coords);
582
583 nr_g = U_DOFS * NODES_PER_EL;
584 nc_g = P_DOFS * NODES_PER_EL;
585
586 for (i = 0; i < nr_g; i++) {
587 for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
588 }
589 }
590
FormStabilisationOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])591 static void FormStabilisationOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
592 {
593 PetscInt ngp;
594 PetscScalar gp_xi[GAUSS_POINTS][2];
595 PetscScalar gp_weight[GAUSS_POINTS];
596 PetscInt p, i, j;
597 PetscScalar Ni_p[NODES_PER_EL];
598 PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
599 PetscScalar J_p, fac, eta_avg;
600
601 /* define quadrature rule */
602 ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
603
604 /* evaluate integral */
605 for (p = 0; p < ngp; p++) {
606 ConstructQ12D_Ni(gp_xi[p], Ni_p);
607 ConstructQ12D_GNi(gp_xi[p], GNi_p);
608 ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
609 fac = gp_weight[p] * J_p;
610
611 for (i = 0; i < NODES_PER_EL; i++) {
612 for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = Ke[NODES_PER_EL * i + j] - fac * (Ni_p[i] * Ni_p[j] - 0.0625);
613 }
614 }
615
616 /* scale */
617 eta_avg = 0.0;
618 for (p = 0; p < ngp; p++) eta_avg += eta[p];
619 eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
620 fac = 1.0 / eta_avg;
621 for (i = 0; i < NODES_PER_EL; i++) {
622 for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
623 }
624 }
625
FormScaledMassMatrixOperatorQ1(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])626 static void FormScaledMassMatrixOperatorQ1(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
627 {
628 PetscInt ngp;
629 PetscScalar gp_xi[GAUSS_POINTS][2];
630 PetscScalar gp_weight[GAUSS_POINTS];
631 PetscInt p, i, j;
632 PetscScalar Ni_p[NODES_PER_EL];
633 PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
634 PetscScalar J_p, fac, eta_avg;
635
636 /* define quadrature rule */
637 ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
638
639 /* evaluate integral */
640 for (p = 0; p < ngp; p++) {
641 ConstructQ12D_Ni(gp_xi[p], Ni_p);
642 ConstructQ12D_GNi(gp_xi[p], GNi_p);
643 ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
644 fac = gp_weight[p] * J_p;
645
646 for (i = 0; i < NODES_PER_EL; i++) {
647 for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = Ke[NODES_PER_EL * i + j] - fac * Ni_p[i] * Ni_p[j];
648 }
649 }
650
651 /* scale */
652 eta_avg = 0.0;
653 for (p = 0; p < ngp; p++) eta_avg += eta[p];
654 eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
655 fac = 1.0 / eta_avg;
656 for (i = 0; i < NODES_PER_EL; i++) {
657 for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
658 }
659 }
660
FormMomentumRhsQ1(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])661 static void FormMomentumRhsQ1(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
662 {
663 PetscInt ngp;
664 PetscScalar gp_xi[GAUSS_POINTS][2];
665 PetscScalar gp_weight[GAUSS_POINTS];
666 PetscInt p, i;
667 PetscScalar Ni_p[NODES_PER_EL];
668 PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
669 PetscScalar J_p, fac;
670
671 /* define quadrature rule */
672 ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
673
674 /* evaluate integral */
675 for (p = 0; p < ngp; p++) {
676 ConstructQ12D_Ni(gp_xi[p], Ni_p);
677 ConstructQ12D_GNi(gp_xi[p], GNi_p);
678 ConstructQ12D_GNx(GNi_p, GNx_p, coords, &J_p);
679 fac = gp_weight[p] * J_p;
680
681 for (i = 0; i < NODES_PER_EL; i++) {
682 Fe[NSD * i] += fac * Ni_p[i] * fx[p];
683 Fe[NSD * i + 1] += fac * Ni_p[i] * fy[p];
684 }
685 }
686 }
687
GetElementCoords(DMDACoor2d ** _coords,PetscInt ei,PetscInt ej,PetscScalar el_coords[])688 static PetscErrorCode GetElementCoords(DMDACoor2d **_coords, PetscInt ei, PetscInt ej, PetscScalar el_coords[])
689 {
690 PetscFunctionBeginUser;
691 /* get coords for the element */
692 el_coords[NSD * 0] = _coords[ej][ei].x;
693 el_coords[NSD * 0 + 1] = _coords[ej][ei].y;
694 el_coords[NSD * 1] = _coords[ej + 1][ei].x;
695 el_coords[NSD * 1 + 1] = _coords[ej + 1][ei].y;
696 el_coords[NSD * 2] = _coords[ej + 1][ei + 1].x;
697 el_coords[NSD * 2 + 1] = _coords[ej + 1][ei + 1].y;
698 el_coords[NSD * 3] = _coords[ej][ei + 1].x;
699 el_coords[NSD * 3 + 1] = _coords[ej][ei + 1].y;
700 PetscFunctionReturn(PETSC_SUCCESS);
701 }
702
AssembleA_Stokes(Mat A,DM stokes_da,DM properties_da,Vec properties)703 static PetscErrorCode AssembleA_Stokes(Mat A, DM stokes_da, DM properties_da, Vec properties)
704 {
705 DM cda;
706 Vec coords;
707 DMDACoor2d **_coords;
708 MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
709 MatStencil p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
710 PetscInt sex, sey, mx, my;
711 PetscInt ei, ej;
712 PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
713 PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
714 PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
715 PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
716 PetscScalar el_coords[NODES_PER_EL * NSD];
717 Vec local_properties;
718 GaussPointCoefficients **props;
719 PetscScalar *prop_eta;
720
721 PetscFunctionBeginUser;
722 /* setup for coords */
723 PetscCall(DMGetCoordinateDM(stokes_da, &cda));
724 PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
725 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
726
727 /* setup for coefficients */
728 PetscCall(DMCreateLocalVector(properties_da, &local_properties));
729 PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
730 PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
731 PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));
732
733 PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
734 PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
735 for (ej = sey; ej < sey + my; ej++) {
736 for (ei = sex; ei < sex + mx; ei++) {
737 /* get coords for the element */
738 PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
739
740 /* get coefficients for the element */
741 prop_eta = props[ej][ei].eta;
742
743 /* initialise element stiffness matrix */
744 PetscCall(PetscMemzero(Ae, sizeof(Ae)));
745 PetscCall(PetscMemzero(Ge, sizeof(Ge)));
746 PetscCall(PetscMemzero(De, sizeof(De)));
747 PetscCall(PetscMemzero(Ce, sizeof(Ce)));
748
749 /* form element stiffness matrix */
750 FormStressOperatorQ1(Ae, el_coords, prop_eta);
751 FormGradientOperatorQ1(Ge, el_coords);
752 FormDivergenceOperatorQ1(De, el_coords);
753 FormStabilisationOperatorQ1(Ce, el_coords, prop_eta);
754
755 /* insert element matrix into global matrix */
756 PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
757 PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
758 PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
759 PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
760 PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
761 }
762 }
763 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
764 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
765
766 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
767
768 PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
769 PetscCall(VecDestroy(&local_properties));
770 PetscFunctionReturn(PETSC_SUCCESS);
771 }
772
AssembleA_PCStokes(Mat A,DM stokes_da,DM properties_da,Vec properties)773 static PetscErrorCode AssembleA_PCStokes(Mat A, DM stokes_da, DM properties_da, Vec properties)
774 {
775 DM cda;
776 Vec coords;
777 DMDACoor2d **_coords;
778 MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
779 MatStencil p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
780 PetscInt sex, sey, mx, my;
781 PetscInt ei, ej;
782 PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
783 PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
784 PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
785 PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
786 PetscScalar el_coords[NODES_PER_EL * NSD];
787 Vec local_properties;
788 GaussPointCoefficients **props;
789 PetscScalar *prop_eta;
790
791 PetscFunctionBeginUser;
792 /* setup for coords */
793 PetscCall(DMGetCoordinateDM(stokes_da, &cda));
794 PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
795 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
796
797 /* setup for coefficients */
798 PetscCall(DMCreateLocalVector(properties_da, &local_properties));
799 PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
800 PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
801 PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));
802
803 PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
804 PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
805 for (ej = sey; ej < sey + my; ej++) {
806 for (ei = sex; ei < sex + mx; ei++) {
807 /* get coords for the element */
808 PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
809
810 /* get coefficients for the element */
811 prop_eta = props[ej][ei].eta;
812
813 /* initialise element stiffness matrix */
814 PetscCall(PetscMemzero(Ae, sizeof(Ae)));
815 PetscCall(PetscMemzero(Ge, sizeof(Ge)));
816 PetscCall(PetscMemzero(De, sizeof(De)));
817 PetscCall(PetscMemzero(Ce, sizeof(Ce)));
818
819 /* form element stiffness matrix */
820 FormStressOperatorQ1(Ae, el_coords, prop_eta);
821 FormGradientOperatorQ1(Ge, el_coords);
822 FormScaledMassMatrixOperatorQ1(Ce, el_coords, prop_eta);
823
824 /* insert element matrix into global matrix */
825 PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
826 PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
827 PetscCall(MatSetValuesStencil(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
828 PetscCall(MatSetValuesStencil(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
829 }
830 }
831 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
832 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
833
834 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
835
836 PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
837 PetscCall(VecDestroy(&local_properties));
838 PetscFunctionReturn(PETSC_SUCCESS);
839 }
840
DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF ** fields_F,MatStencil u_eqn[],MatStencil p_eqn[],PetscScalar Fe_u[],PetscScalar Fe_p[])841 static PetscErrorCode DMDASetValuesLocalStencil_ADD_VALUES(StokesDOF **fields_F, MatStencil u_eqn[], MatStencil p_eqn[], PetscScalar Fe_u[], PetscScalar Fe_p[])
842 {
843 PetscInt n;
844
845 PetscFunctionBeginUser;
846 for (n = 0; n < 4; n++) {
847 fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].u_dof = fields_F[u_eqn[2 * n].j][u_eqn[2 * n].i].u_dof + Fe_u[2 * n];
848 fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].v_dof = fields_F[u_eqn[2 * n + 1].j][u_eqn[2 * n + 1].i].v_dof + Fe_u[2 * n + 1];
849 fields_F[p_eqn[n].j][p_eqn[n].i].p_dof = fields_F[p_eqn[n].j][p_eqn[n].i].p_dof + Fe_p[n];
850 }
851 PetscFunctionReturn(PETSC_SUCCESS);
852 }
853
AssembleF_Stokes(Vec F,DM stokes_da,DM properties_da,Vec properties)854 static PetscErrorCode AssembleF_Stokes(Vec F, DM stokes_da, DM properties_da, Vec properties)
855 {
856 DM cda;
857 Vec coords;
858 DMDACoor2d **_coords;
859 MatStencil u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
860 MatStencil p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
861 PetscInt sex, sey, mx, my;
862 PetscInt ei, ej;
863 PetscScalar Fe[NODES_PER_EL * U_DOFS];
864 PetscScalar He[NODES_PER_EL * P_DOFS];
865 PetscScalar el_coords[NODES_PER_EL * NSD];
866 Vec local_properties;
867 GaussPointCoefficients **props;
868 PetscScalar *prop_fx, *prop_fy;
869 Vec local_F;
870 StokesDOF **ff;
871
872 PetscFunctionBeginUser;
873 /* setup for coords */
874 PetscCall(DMGetCoordinateDM(stokes_da, &cda));
875 PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
876 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
877
878 /* setup for coefficients */
879 PetscCall(DMGetLocalVector(properties_da, &local_properties));
880 PetscCall(DMGlobalToLocalBegin(properties_da, properties, INSERT_VALUES, local_properties));
881 PetscCall(DMGlobalToLocalEnd(properties_da, properties, INSERT_VALUES, local_properties));
882 PetscCall(DMDAVecGetArray(properties_da, local_properties, &props));
883
884 /* get access to the vector */
885 PetscCall(DMGetLocalVector(stokes_da, &local_F));
886 PetscCall(VecZeroEntries(local_F));
887 PetscCall(DMDAVecGetArray(stokes_da, local_F, &ff));
888
889 PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
890 PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
891 for (ej = sey; ej < sey + my; ej++) {
892 for (ei = sex; ei < sex + mx; ei++) {
893 /* get coords for the element */
894 PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
895
896 /* get coefficients for the element */
897 prop_fx = props[ej][ei].fx;
898 prop_fy = props[ej][ei].fy;
899
900 /* initialise element stiffness matrix */
901 PetscCall(PetscMemzero(Fe, sizeof(Fe)));
902 PetscCall(PetscMemzero(He, sizeof(He)));
903
904 /* form element stiffness matrix */
905 FormMomentumRhsQ1(Fe, el_coords, prop_fx, prop_fy);
906
907 /* insert element matrix into global matrix */
908 PetscCall(DMDAGetElementEqnums_up(u_eqn, p_eqn, ei, ej));
909
910 PetscCall(DMDASetValuesLocalStencil_ADD_VALUES(ff, u_eqn, p_eqn, Fe, He));
911 }
912 }
913
914 PetscCall(DMDAVecRestoreArray(stokes_da, local_F, &ff));
915 PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
916 PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
917 PetscCall(DMRestoreLocalVector(stokes_da, &local_F));
918
919 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
920
921 PetscCall(DMDAVecRestoreArray(properties_da, local_properties, &props));
922 PetscCall(DMRestoreLocalVector(properties_da, &local_properties));
923 PetscFunctionReturn(PETSC_SUCCESS);
924 }
925
DMDACreateSolCx(PetscReal eta0,PetscReal eta1,PetscReal xc,PetscInt nz,PetscInt mx,PetscInt my,DM * _da,Vec * _X)926 static PetscErrorCode DMDACreateSolCx(PetscReal eta0, PetscReal eta1, PetscReal xc, PetscInt nz, PetscInt mx, PetscInt my, DM *_da, Vec *_X)
927 {
928 DM da, cda;
929 Vec X;
930 StokesDOF **_stokes;
931 Vec coords;
932 DMDACoor2d **_coords;
933 PetscInt si, sj, ei, ej, i, j;
934
935 PetscFunctionBeginUser;
936 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, 3, 1, NULL, NULL, &da));
937 PetscCall(DMSetFromOptions(da));
938 PetscCall(DMSetUp(da));
939 PetscCall(DMDASetFieldName(da, 0, "anlytic_Vx"));
940 PetscCall(DMDASetFieldName(da, 1, "anlytic_Vy"));
941 PetscCall(DMDASetFieldName(da, 2, "analytic_P"));
942
943 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0., 0.));
944
945 PetscCall(DMGetCoordinatesLocal(da, &coords));
946 PetscCall(DMGetCoordinateDM(da, &cda));
947 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
948
949 PetscCall(DMCreateGlobalVector(da, &X));
950 PetscCall(DMDAVecGetArray(da, X, &_stokes));
951
952 PetscCall(DMDAGetCorners(da, &si, &sj, 0, &ei, &ej, 0));
953 for (j = sj; j < sj + ej; j++) {
954 for (i = si; i < si + ei; i++) {
955 PetscReal pos[2], pressure, vel[2], total_stress[3], strain_rate[3];
956
957 pos[0] = PetscRealPart(_coords[j][i].x);
958 pos[1] = PetscRealPart(_coords[j][i].y);
959
960 evaluate_solCx(pos, eta0, eta1, xc, nz, vel, &pressure, total_stress, strain_rate);
961
962 _stokes[j][i].u_dof = vel[0];
963 _stokes[j][i].v_dof = vel[1];
964 _stokes[j][i].p_dof = pressure;
965 }
966 }
967 PetscCall(DMDAVecRestoreArray(da, X, &_stokes));
968 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
969
970 *_da = da;
971 *_X = X;
972 PetscFunctionReturn(PETSC_SUCCESS);
973 }
974
StokesDAGetNodalFields(StokesDOF ** fields,PetscInt ei,PetscInt ej,StokesDOF nodal_fields[])975 static PetscErrorCode StokesDAGetNodalFields(StokesDOF **fields, PetscInt ei, PetscInt ej, StokesDOF nodal_fields[])
976 {
977 PetscFunctionBeginUser;
978 /* get the nodal fields */
979 nodal_fields[0].u_dof = fields[ej][ei].u_dof;
980 nodal_fields[0].v_dof = fields[ej][ei].v_dof;
981 nodal_fields[0].p_dof = fields[ej][ei].p_dof;
982 nodal_fields[1].u_dof = fields[ej + 1][ei].u_dof;
983 nodal_fields[1].v_dof = fields[ej + 1][ei].v_dof;
984 nodal_fields[1].p_dof = fields[ej + 1][ei].p_dof;
985 nodal_fields[2].u_dof = fields[ej + 1][ei + 1].u_dof;
986 nodal_fields[2].v_dof = fields[ej + 1][ei + 1].v_dof;
987 nodal_fields[2].p_dof = fields[ej + 1][ei + 1].p_dof;
988 nodal_fields[3].u_dof = fields[ej][ei + 1].u_dof;
989 nodal_fields[3].v_dof = fields[ej][ei + 1].v_dof;
990 nodal_fields[3].p_dof = fields[ej][ei + 1].p_dof;
991 PetscFunctionReturn(PETSC_SUCCESS);
992 }
993
DMDAIntegrateErrors(DM stokes_da,Vec X,Vec X_analytic)994 static PetscErrorCode DMDAIntegrateErrors(DM stokes_da, Vec X, Vec X_analytic)
995 {
996 DM cda;
997 Vec coords, X_analytic_local, X_local;
998 DMDACoor2d **_coords;
999 PetscInt sex, sey, mx, my;
1000 PetscInt ei, ej;
1001 PetscScalar el_coords[NODES_PER_EL * NSD];
1002 StokesDOF **stokes_analytic, **stokes;
1003 StokesDOF stokes_analytic_e[4], stokes_e[4];
1004
1005 PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
1006 PetscScalar Ni_p[NODES_PER_EL];
1007 PetscInt ngp;
1008 PetscScalar gp_xi[GAUSS_POINTS][2];
1009 PetscScalar gp_weight[GAUSS_POINTS];
1010 PetscInt p, i;
1011 PetscScalar J_p, fac;
1012 PetscScalar h, p_e_L2, u_e_L2, u_e_H1, p_L2, u_L2, u_H1, tp_L2, tu_L2, tu_H1;
1013 PetscInt M;
1014 PetscReal xymin[2], xymax[2];
1015
1016 PetscFunctionBeginUser;
1017 /* define quadrature rule */
1018 ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
1019
1020 /* setup for coords */
1021 PetscCall(DMGetCoordinateDM(stokes_da, &cda));
1022 PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
1023 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1024
1025 /* setup for analytic */
1026 PetscCall(DMCreateLocalVector(stokes_da, &X_analytic_local));
1027 PetscCall(DMGlobalToLocalBegin(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
1028 PetscCall(DMGlobalToLocalEnd(stokes_da, X_analytic, INSERT_VALUES, X_analytic_local));
1029 PetscCall(DMDAVecGetArray(stokes_da, X_analytic_local, &stokes_analytic));
1030
1031 /* setup for solution */
1032 PetscCall(DMCreateLocalVector(stokes_da, &X_local));
1033 PetscCall(DMGlobalToLocalBegin(stokes_da, X, INSERT_VALUES, X_local));
1034 PetscCall(DMGlobalToLocalEnd(stokes_da, X, INSERT_VALUES, X_local));
1035 PetscCall(DMDAVecGetArray(stokes_da, X_local, &stokes));
1036
1037 PetscCall(DMDAGetInfo(stokes_da, 0, &M, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1038 PetscCall(DMGetBoundingBox(stokes_da, xymin, xymax));
1039
1040 h = (xymax[0] - xymin[0]) / ((PetscReal)M);
1041
1042 tp_L2 = tu_L2 = tu_H1 = 0.0;
1043
1044 PetscCall(DMDAGetElementsCorners(stokes_da, &sex, &sey, NULL));
1045 PetscCall(DMDAGetElementsSizes(stokes_da, &mx, &my, NULL));
1046 for (ej = sey; ej < sey + my; ej++) {
1047 for (ei = sex; ei < sex + mx; ei++) {
1048 /* get coords for the element */
1049 PetscCall(GetElementCoords(_coords, ei, ej, el_coords));
1050 PetscCall(StokesDAGetNodalFields(stokes, ei, ej, stokes_e));
1051 PetscCall(StokesDAGetNodalFields(stokes_analytic, ei, ej, stokes_analytic_e));
1052
1053 /* evaluate integral */
1054 p_e_L2 = 0.0;
1055 u_e_L2 = 0.0;
1056 u_e_H1 = 0.0;
1057 for (p = 0; p < ngp; p++) {
1058 ConstructQ12D_Ni(gp_xi[p], Ni_p);
1059 ConstructQ12D_GNi(gp_xi[p], GNi_p);
1060 ConstructQ12D_GNx(GNi_p, GNx_p, el_coords, &J_p);
1061 fac = gp_weight[p] * J_p;
1062
1063 for (i = 0; i < NODES_PER_EL; i++) {
1064 PetscScalar u_error, v_error;
1065
1066 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);
1067
1068 u_error = stokes_e[i].u_dof - stokes_analytic_e[i].u_dof;
1069 v_error = stokes_e[i].v_dof - stokes_analytic_e[i].v_dof;
1070 u_e_L2 += fac * Ni_p[i] * (u_error * u_error + v_error * v_error);
1071
1072 u_e_H1 = u_e_H1 + fac * (GNx_p[0][i] * u_error * GNx_p[0][i] * u_error /* du/dx */
1073 + GNx_p[1][i] * u_error * GNx_p[1][i] * u_error /* du/dy */
1074 + GNx_p[0][i] * v_error * GNx_p[0][i] * v_error /* dv/dx */
1075 + GNx_p[1][i] * v_error * GNx_p[1][i] * v_error); /* dv/dy */
1076 }
1077 }
1078
1079 tp_L2 += p_e_L2;
1080 tu_L2 += u_e_L2;
1081 tu_H1 += u_e_H1;
1082 }
1083 }
1084 PetscCallMPI(MPIU_Allreduce(&tp_L2, &p_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1085 PetscCallMPI(MPIU_Allreduce(&tu_L2, &u_L2, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1086 PetscCallMPI(MPIU_Allreduce(&tu_H1, &u_H1, 1, MPIU_SCALAR, MPIU_SUM, PETSC_COMM_WORLD));
1087 p_L2 = PetscSqrtScalar(p_L2);
1088 u_L2 = PetscSqrtScalar(u_L2);
1089 u_H1 = PetscSqrtScalar(u_H1);
1090
1091 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%1.4e %1.4e %1.4e %1.4e\n", (double)PetscRealPart(h), (double)PetscRealPart(p_L2), (double)PetscRealPart(u_L2), (double)PetscRealPart(u_H1)));
1092
1093 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1094
1095 PetscCall(DMDAVecRestoreArray(stokes_da, X_analytic_local, &stokes_analytic));
1096 PetscCall(VecDestroy(&X_analytic_local));
1097 PetscCall(DMDAVecRestoreArray(stokes_da, X_local, &stokes));
1098 PetscCall(VecDestroy(&X_local));
1099 PetscFunctionReturn(PETSC_SUCCESS);
1100 }
1101
solve_stokes_2d_coupled(PetscInt mx,PetscInt my)1102 static PetscErrorCode solve_stokes_2d_coupled(PetscInt mx, PetscInt my)
1103 {
1104 DM da_Stokes, da_prop;
1105 PetscInt u_dof, p_dof, dof, stencil_width;
1106 Mat A, B;
1107 DM prop_cda, vel_cda;
1108 Vec prop_coords, vel_coords;
1109 PetscInt si, sj, nx, ny, i, j, p;
1110 Vec f, X;
1111 PetscInt prop_dof, prop_stencil_width;
1112 Vec properties, l_properties;
1113 PetscReal dx, dy;
1114 PetscInt M, N;
1115 DMDACoor2d **_prop_coords, **_vel_coords;
1116 GaussPointCoefficients **element_props;
1117 PetscInt its;
1118 KSP ksp_S;
1119 PetscInt coefficient_structure = 0;
1120 PetscInt cpu_x, cpu_y, *lx = NULL, *ly = NULL;
1121 PetscBool use_gp_coords = PETSC_FALSE, set, output_gnuplot = PETSC_FALSE, glvis = PETSC_FALSE, change = PETSC_FALSE;
1122 char filename[PETSC_MAX_PATH_LEN];
1123
1124 PetscFunctionBeginUser;
1125 PetscCall(PetscOptionsGetBool(NULL, NULL, "-gnuplot", &output_gnuplot, NULL));
1126 PetscCall(PetscOptionsGetBool(NULL, NULL, "-glvis", &glvis, NULL));
1127 PetscCall(PetscOptionsGetBool(NULL, NULL, "-change", &change, NULL));
1128
1129 /* Generate the da for velocity and pressure */
1130 /*
1131 We use Q1 elements for the temperature.
1132 FEM has a 9-point stencil (BOX) or connectivity pattern
1133 Num nodes in each direction is mx+1, my+1
1134 */
1135 u_dof = U_DOFS; /* Vx, Vy - velocities */
1136 p_dof = P_DOFS; /* p - pressure */
1137 dof = u_dof + p_dof;
1138 stencil_width = 1;
1139 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, &da_Stokes));
1140
1141 PetscCall(DMSetMatType(da_Stokes, MATAIJ));
1142 PetscCall(DMSetFromOptions(da_Stokes));
1143 PetscCall(DMSetUp(da_Stokes));
1144 PetscCall(DMDASetFieldName(da_Stokes, 0, "Vx"));
1145 PetscCall(DMDASetFieldName(da_Stokes, 1, "Vy"));
1146 PetscCall(DMDASetFieldName(da_Stokes, 2, "P"));
1147
1148 /* unit box [0,1] x [0,1] */
1149 PetscCall(DMDASetUniformCoordinates(da_Stokes, 0.0, 1.0, 0.0, 1.0, 0., 0.));
1150
1151 /* Generate element properties, we will assume all material properties are constant over the element */
1152 /* !!! IN PARALLEL WE MUST MAKE SURE THE TWO DMDA's ALIGN !!! */
1153 PetscCall(DMDAGetInfo(da_Stokes, 0, 0, 0, 0, &cpu_x, &cpu_y, 0, 0, 0, 0, 0, 0, 0));
1154 PetscCall(DMDAGetElementOwnershipRanges2d(da_Stokes, &lx, &ly));
1155
1156 prop_dof = (PetscInt)(sizeof(GaussPointCoefficients) / sizeof(PetscScalar)); /* gauss point setup */
1157 prop_stencil_width = 0;
1158 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx, my, cpu_x, cpu_y, prop_dof, prop_stencil_width, lx, ly, &da_prop));
1159 PetscCall(DMSetFromOptions(da_prop));
1160 PetscCall(DMSetUp(da_prop));
1161 PetscCall(PetscFree(lx));
1162 PetscCall(PetscFree(ly));
1163
1164 /* define centroid positions */
1165 PetscCall(DMDAGetInfo(da_prop, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0));
1166 dx = 1.0 / (PetscReal)M;
1167 dy = 1.0 / (PetscReal)N;
1168
1169 PetscCall(DMDASetUniformCoordinates(da_prop, 0.0 + 0.5 * dx, 1.0 - 0.5 * dx, 0.0 + 0.5 * dy, 1.0 - 0.5 * dy, 0., 0));
1170
1171 /* define coefficients */
1172 PetscCall(PetscOptionsGetInt(NULL, NULL, "-c_str", &coefficient_structure, NULL));
1173
1174 PetscCall(DMCreateGlobalVector(da_prop, &properties));
1175 PetscCall(DMCreateLocalVector(da_prop, &l_properties));
1176 PetscCall(DMDAVecGetArray(da_prop, l_properties, &element_props));
1177
1178 PetscCall(DMGetCoordinateDM(da_prop, &prop_cda));
1179 PetscCall(DMGetCoordinatesLocal(da_prop, &prop_coords));
1180 PetscCall(DMDAVecGetArray(prop_cda, prop_coords, &_prop_coords));
1181
1182 PetscCall(DMDAGetGhostCorners(prop_cda, &si, &sj, 0, &nx, &ny, 0));
1183
1184 PetscCall(DMGetCoordinateDM(da_Stokes, &vel_cda));
1185 PetscCall(DMGetCoordinatesLocal(da_Stokes, &vel_coords));
1186 PetscCall(DMDAVecGetArray(vel_cda, vel_coords, &_vel_coords));
1187
1188 /* interpolate the coordinates */
1189 for (j = sj; j < sj + ny; j++) {
1190 for (i = si; i < si + nx; i++) {
1191 PetscInt ngp;
1192 PetscScalar gp_xi[GAUSS_POINTS][2], gp_weight[GAUSS_POINTS];
1193 PetscScalar el_coords[8];
1194
1195 PetscCall(GetElementCoords(_vel_coords, i, j, el_coords));
1196 ConstructGaussQuadrature(&ngp, gp_xi, gp_weight);
1197
1198 for (p = 0; p < GAUSS_POINTS; p++) {
1199 PetscScalar gp_x, gp_y;
1200 PetscInt n;
1201 PetscScalar xi_p[2], Ni_p[4];
1202
1203 xi_p[0] = gp_xi[p][0];
1204 xi_p[1] = gp_xi[p][1];
1205 ConstructQ12D_Ni(xi_p, Ni_p);
1206
1207 gp_x = 0.0;
1208 gp_y = 0.0;
1209 for (n = 0; n < NODES_PER_EL; n++) {
1210 gp_x = gp_x + Ni_p[n] * el_coords[2 * n];
1211 gp_y = gp_y + Ni_p[n] * el_coords[2 * n + 1];
1212 }
1213 element_props[j][i].gp_coords[2 * p] = gp_x;
1214 element_props[j][i].gp_coords[2 * p + 1] = gp_y;
1215 }
1216 }
1217 }
1218
1219 /* define the coefficients */
1220 PetscCall(PetscOptionsGetBool(NULL, NULL, "-use_gp_coords", &use_gp_coords, NULL));
1221
1222 for (j = sj; j < sj + ny; j++) {
1223 for (i = si; i < si + nx; i++) {
1224 PetscReal centroid_x = PetscRealPart(_prop_coords[j][i].x); /* centroids of cell */
1225 PetscReal centroid_y = PetscRealPart(_prop_coords[j][i].y);
1226 PetscReal coord_x, coord_y;
1227
1228 if (coefficient_structure == 0) {
1229 PetscReal opts_eta0, opts_eta1, opts_xc;
1230 PetscInt opts_nz;
1231
1232 opts_eta0 = 1.0;
1233 opts_eta1 = 1.0;
1234 opts_xc = 0.5;
1235 opts_nz = 1;
1236
1237 PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL));
1238 PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL));
1239 PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL));
1240 PetscCall(PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL));
1241
1242 for (p = 0; p < GAUSS_POINTS; p++) {
1243 coord_x = centroid_x;
1244 coord_y = centroid_y;
1245 if (use_gp_coords) {
1246 coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1247 coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1248 }
1249
1250 element_props[j][i].eta[p] = opts_eta0;
1251 if (coord_x > opts_xc) element_props[j][i].eta[p] = opts_eta1;
1252
1253 element_props[j][i].fx[p] = 0.0;
1254 element_props[j][i].fy[p] = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1255 }
1256 } else if (coefficient_structure == 1) { /* square sinker */
1257 PetscReal opts_eta0, opts_eta1, opts_dx, opts_dy;
1258
1259 opts_eta0 = 1.0;
1260 opts_eta1 = 1.0;
1261 opts_dx = 0.50;
1262 opts_dy = 0.50;
1263
1264 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1265 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1266 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL));
1267 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL));
1268
1269 for (p = 0; p < GAUSS_POINTS; p++) {
1270 coord_x = centroid_x;
1271 coord_y = centroid_y;
1272 if (use_gp_coords) {
1273 coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1274 coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1275 }
1276
1277 element_props[j][i].eta[p] = opts_eta0;
1278 element_props[j][i].fx[p] = 0.0;
1279 element_props[j][i].fy[p] = 0.0;
1280
1281 if ((coord_x > -0.5 * opts_dx + 0.5) && (coord_x < 0.5 * opts_dx + 0.5)) {
1282 if ((coord_y > -0.5 * opts_dy + 0.5) && (coord_y < 0.5 * opts_dy + 0.5)) {
1283 element_props[j][i].eta[p] = opts_eta1;
1284 element_props[j][i].fx[p] = 0.0;
1285 element_props[j][i].fy[p] = -1.0;
1286 }
1287 }
1288 }
1289 } else if (coefficient_structure == 2) { /* circular sinker */
1290 PetscReal opts_eta0, opts_eta1, opts_r, radius2;
1291
1292 opts_eta0 = 1.0;
1293 opts_eta1 = 1.0;
1294 opts_r = 0.25;
1295
1296 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1297 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1298 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL));
1299
1300 for (p = 0; p < GAUSS_POINTS; p++) {
1301 coord_x = centroid_x;
1302 coord_y = centroid_y;
1303 if (use_gp_coords) {
1304 coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1305 coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1306 }
1307
1308 element_props[j][i].eta[p] = opts_eta0;
1309 element_props[j][i].fx[p] = 0.0;
1310 element_props[j][i].fy[p] = 0.0;
1311
1312 radius2 = (coord_x - 0.5) * (coord_x - 0.5) + (coord_y - 0.5) * (coord_y - 0.5);
1313 if (radius2 < opts_r * opts_r) {
1314 element_props[j][i].eta[p] = opts_eta1;
1315 element_props[j][i].fx[p] = 0.0;
1316 element_props[j][i].fy[p] = -1.0;
1317 }
1318 }
1319 } else if (coefficient_structure == 3) { /* circular and rectangular inclusion */
1320 PetscReal opts_eta0, opts_eta1, opts_r, opts_dx, opts_dy, opts_c0x, opts_c0y, opts_s0x, opts_s0y, opts_phi, radius2;
1321
1322 opts_eta0 = 1.0;
1323 opts_eta1 = 1.0;
1324 opts_r = 0.25;
1325 opts_c0x = 0.35; /* circle center */
1326 opts_c0y = 0.35;
1327 opts_s0x = 0.7; /* square center */
1328 opts_s0y = 0.7;
1329 opts_dx = 0.25;
1330 opts_dy = 0.25;
1331 opts_phi = 25;
1332
1333 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta0", &opts_eta0, NULL));
1334 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_eta1", &opts_eta1, NULL));
1335 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_r", &opts_r, NULL));
1336 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_c0x", &opts_c0x, NULL));
1337 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_c0y", &opts_c0y, NULL));
1338 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_s0x", &opts_s0x, NULL));
1339 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_s0y", &opts_s0y, NULL));
1340 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dx", &opts_dx, NULL));
1341 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_dy", &opts_dy, NULL));
1342 PetscCall(PetscOptionsGetReal(NULL, NULL, "-sinker_phi", &opts_phi, NULL));
1343 opts_phi *= PETSC_PI / 180;
1344
1345 for (p = 0; p < GAUSS_POINTS; p++) {
1346 coord_x = centroid_x;
1347 coord_y = centroid_y;
1348 if (use_gp_coords) {
1349 coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1350 coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1351 }
1352
1353 element_props[j][i].eta[p] = opts_eta0;
1354 element_props[j][i].fx[p] = 0.0;
1355 element_props[j][i].fy[p] = -0.2;
1356
1357 radius2 = PetscSqr(coord_x - opts_c0x) + PetscSqr(coord_y - opts_c0y);
1358 if (radius2 < opts_r * opts_r || (PetscAbs(+(coord_x - opts_s0x) * PetscCosReal(opts_phi) + (coord_y - opts_s0y) * PetscSinReal(opts_phi)) < opts_dx / 2 && PetscAbs(-(coord_x - opts_s0x) * PetscSinReal(opts_phi) + (coord_y - opts_s0y) * PetscCosReal(opts_phi)) < opts_dy / 2)) {
1359 element_props[j][i].eta[p] = opts_eta1;
1360 element_props[j][i].fx[p] = 0.0;
1361 element_props[j][i].fy[p] = -1.0;
1362 }
1363 }
1364 } else if (coefficient_structure == 4) { /* subdomain jump */
1365 PetscReal opts_mag, opts_eta0;
1366 PetscInt opts_nz, px, py;
1367 PetscBool jump;
1368
1369 opts_mag = 1.0;
1370 opts_eta0 = 1.0;
1371 opts_nz = 1;
1372
1373 PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_eta0", &opts_eta0, NULL));
1374 PetscCall(PetscOptionsGetReal(NULL, NULL, "-jump_magnitude", &opts_mag, NULL));
1375 PetscCall(PetscOptionsGetInt(NULL, NULL, "-jump_nz", &opts_nz, NULL));
1376 PetscCall(DMDAGetInfo(da_Stokes, NULL, NULL, NULL, NULL, &px, &py, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
1377 if (px % 2) {
1378 jump = (PetscBool)(PetscGlobalRank % 2);
1379 } else {
1380 jump = (PetscBool)((PetscGlobalRank / px) % 2 ? PetscGlobalRank % 2 : !(PetscGlobalRank % 2));
1381 }
1382 for (p = 0; p < GAUSS_POINTS; p++) {
1383 coord_x = centroid_x;
1384 coord_y = centroid_y;
1385 if (use_gp_coords) {
1386 coord_x = PetscRealPart(element_props[j][i].gp_coords[2 * p]);
1387 coord_y = PetscRealPart(element_props[j][i].gp_coords[2 * p + 1]);
1388 }
1389
1390 element_props[j][i].eta[p] = jump ? PetscPowReal(10.0, opts_mag) : opts_eta0;
1391 element_props[j][i].fx[p] = 0.0;
1392 element_props[j][i].fy[p] = PetscSinReal(opts_nz * PETSC_PI * coord_y) * PetscCosReal(1.0 * PETSC_PI * coord_x);
1393 }
1394 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_USER, "Unknown coefficient_structure");
1395 }
1396 }
1397 PetscCall(DMDAVecRestoreArray(prop_cda, prop_coords, &_prop_coords));
1398
1399 PetscCall(DMDAVecRestoreArray(vel_cda, vel_coords, &_vel_coords));
1400
1401 PetscCall(DMDAVecRestoreArray(da_prop, l_properties, &element_props));
1402 PetscCall(DMLocalToGlobalBegin(da_prop, l_properties, ADD_VALUES, properties));
1403 PetscCall(DMLocalToGlobalEnd(da_prop, l_properties, ADD_VALUES, properties));
1404
1405 if (output_gnuplot) {
1406 PetscCall(DMDACoordViewGnuplot2d(da_Stokes, "mesh"));
1407 PetscCall(DMDAViewCoefficientsGnuplot2d(da_prop, properties, "Coefficients for Stokes eqn.", "properties"));
1408 }
1409
1410 if (glvis) {
1411 Vec glv_prop, etaf;
1412 PetscViewer view;
1413 PetscInt dim;
1414 const char *fec = {"FiniteElementCollection: L2_2D_P1"};
1415
1416 PetscCall(DMGetDimension(da_Stokes, &dim));
1417 PetscCall(VecCreateSeq(PETSC_COMM_SELF, GAUSS_POINTS * mx * mx, &etaf));
1418 PetscCall(PetscObjectSetName((PetscObject)etaf, "viscosity"));
1419 PetscCall(PetscViewerGLVisOpen(PETSC_COMM_WORLD, PETSC_VIEWER_GLVIS_SOCKET, NULL, PETSC_DECIDE, &view));
1420 PetscCall(PetscViewerGLVisSetFields(view, 1, &fec, &dim, glvis_extract_eta, (PetscObject *)&etaf, da_prop, NULL));
1421 PetscCall(DMGetLocalVector(da_prop, &glv_prop));
1422 PetscCall(DMGlobalToLocalBegin(da_prop, properties, INSERT_VALUES, glv_prop));
1423 PetscCall(DMGlobalToLocalEnd(da_prop, properties, INSERT_VALUES, glv_prop));
1424 PetscCall(VecSetDM(etaf, da_Stokes));
1425 PetscCall(VecView(glv_prop, view));
1426 PetscCall(PetscViewerDestroy(&view));
1427 PetscCall(DMRestoreLocalVector(da_prop, &glv_prop));
1428 PetscCall(VecDestroy(&etaf));
1429 }
1430
1431 /* Generate a matrix with the correct non-zero pattern of type AIJ. This will work in parallel and serial */
1432 PetscCall(DMCreateMatrix(da_Stokes, &A));
1433 PetscCall(DMCreateMatrix(da_Stokes, &B));
1434 PetscCall(DMCreateGlobalVector(da_Stokes, &f));
1435 PetscCall(DMCreateGlobalVector(da_Stokes, &X));
1436
1437 /* assemble A11 */
1438 PetscCall(MatZeroEntries(A));
1439 PetscCall(MatZeroEntries(B));
1440 PetscCall(VecZeroEntries(f));
1441
1442 PetscCall(AssembleA_Stokes(A, da_Stokes, da_prop, properties));
1443 PetscCall(AssembleA_PCStokes(B, da_Stokes, da_prop, properties));
1444 /* build force vector */
1445 PetscCall(AssembleF_Stokes(f, da_Stokes, da_prop, properties));
1446
1447 PetscCall(DMDABCApplyFreeSlip(da_Stokes, A, f));
1448 PetscCall(DMDABCApplyFreeSlip(da_Stokes, B, NULL));
1449
1450 /* SOLVE */
1451 PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp_S));
1452 PetscCall(KSPSetOptionsPrefix(ksp_S, "stokes_"));
1453 PetscCall(KSPSetDM(ksp_S, da_Stokes));
1454 PetscCall(KSPSetDMActive(ksp_S, KSP_DMACTIVE_ALL, PETSC_FALSE));
1455 PetscCall(KSPSetOperators(ksp_S, A, B));
1456 PetscCall(KSPSetFromOptions(ksp_S));
1457 {
1458 PC pc;
1459 const PetscInt ufields[] = {0, 1}, pfields[1] = {2};
1460
1461 PetscCall(KSPGetPC(ksp_S, &pc));
1462 PetscCall(PCFieldSplitSetBlockSize(pc, 3));
1463 PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
1464 PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
1465 }
1466
1467 {
1468 PC pc;
1469 PetscBool same = PETSC_FALSE;
1470 PetscCall(KSPGetPC(ksp_S, &pc));
1471 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &same));
1472 if (same) {
1473 PetscBool usedivmat = PETSC_FALSE;
1474 PetscCall(KSPSetOperators(ksp_S, A, A));
1475
1476 PetscCall(PetscOptionsGetBool(NULL, NULL, "-stokes_pc_bddc_use_divergence", &usedivmat, NULL));
1477 if (usedivmat) {
1478 IS *fields, vel;
1479 PetscInt i, nf;
1480
1481 PetscCall(DMCreateFieldDecomposition(da_Stokes, &nf, NULL, &fields, NULL));
1482 PetscCall(ISConcatenate(PETSC_COMM_WORLD, 2, fields, &vel));
1483 PetscCall(MatZeroRowsIS(B, fields[2], 1.0, NULL, NULL)); /* we put 1.0 on the diagonal to pick the pressure average too */
1484 PetscCall(MatTranspose(B, MAT_INPLACE_MATRIX, &B));
1485 PetscCall(MatZeroRowsIS(B, vel, 0.0, NULL, NULL));
1486 PetscCall(ISDestroy(&vel));
1487 PetscCall(PCBDDCSetDivergenceMat(pc, B, PETSC_FALSE, NULL));
1488 for (i = 0; i < nf; i++) PetscCall(ISDestroy(&fields[i]));
1489 PetscCall(PetscFree(fields));
1490 }
1491 }
1492 }
1493
1494 {
1495 PC pc_S;
1496 KSP *sub_ksp, ksp_U;
1497 PetscInt nsplits;
1498 DM da_U;
1499 PetscBool is_pcfs;
1500
1501 PetscCall(KSPSetUp(ksp_S));
1502 PetscCall(KSPGetPC(ksp_S, &pc_S));
1503
1504 is_pcfs = PETSC_FALSE;
1505 PetscCall(PetscObjectTypeCompare((PetscObject)pc_S, PCFIELDSPLIT, &is_pcfs));
1506
1507 if (is_pcfs) {
1508 PetscCall(PCFieldSplitGetSubKSP(pc_S, &nsplits, &sub_ksp));
1509 ksp_U = sub_ksp[0];
1510 PetscCall(PetscFree(sub_ksp));
1511
1512 if (nsplits == 2) {
1513 PetscCall(DMDACreateCompatibleDMDA(da_Stokes, 2, &da_U));
1514
1515 PetscCall(KSPSetDM(ksp_U, da_U));
1516 PetscCall(KSPSetDMActive(ksp_U, KSP_DMACTIVE_ALL, PETSC_FALSE));
1517 PetscCall(DMDestroy(&da_U));
1518 if (change) {
1519 const char opt[] = "-stokes_fieldsplit_u_pc_factor_mat_solver_type mumps";
1520 PetscCall(PetscOptionsInsertString(NULL, opt));
1521 PetscCall(KSPSetFromOptions(ksp_U));
1522 }
1523 PetscCall(KSPSetUp(ksp_U));
1524 }
1525 }
1526 }
1527
1528 PetscCall(KSPSolve(ksp_S, f, X));
1529
1530 PetscCall(PetscOptionsGetString(NULL, NULL, "-o", filename, sizeof(filename), &set));
1531 if (set) {
1532 char *ext;
1533 PetscViewer viewer;
1534 PetscBool flg;
1535 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1536 PetscCall(PetscStrrchr(filename, '.', &ext));
1537 PetscCall(PetscStrcmp("vts", ext, &flg));
1538 if (flg) {
1539 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1540 } else {
1541 PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
1542 }
1543 PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1544 PetscCall(PetscViewerFileSetName(viewer, filename));
1545 PetscCall(VecView(X, viewer));
1546 PetscCall(PetscViewerDestroy(&viewer));
1547 }
1548 if (output_gnuplot) PetscCall(DMDAViewGnuplot2d(da_Stokes, X, "Velocity solution for Stokes eqn.", "X"));
1549
1550 if (glvis) {
1551 PetscViewer view;
1552
1553 PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &view));
1554 PetscCall(PetscViewerSetType(view, PETSCVIEWERGLVIS));
1555 PetscCall(VecView(X, view));
1556 PetscCall(PetscViewerDestroy(&view));
1557 }
1558
1559 PetscCall(KSPGetIterationNumber(ksp_S, &its));
1560
1561 if (coefficient_structure == 0) {
1562 PetscReal opts_eta0, opts_eta1, opts_xc;
1563 PetscInt opts_nz, N;
1564 DM da_Stokes_analytic;
1565 Vec X_analytic;
1566 PetscReal nrm1[3], nrm2[3], nrmI[3];
1567
1568 opts_eta0 = 1.0;
1569 opts_eta1 = 1.0;
1570 opts_xc = 0.5;
1571 opts_nz = 1;
1572
1573 PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta0", &opts_eta0, NULL));
1574 PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_eta1", &opts_eta1, NULL));
1575 PetscCall(PetscOptionsGetReal(NULL, NULL, "-solcx_xc", &opts_xc, NULL));
1576 PetscCall(PetscOptionsGetInt(NULL, NULL, "-solcx_nz", &opts_nz, NULL));
1577
1578 PetscCall(DMDACreateSolCx(opts_eta0, opts_eta1, opts_xc, opts_nz, mx, my, &da_Stokes_analytic, &X_analytic));
1579 if (output_gnuplot) PetscCall(DMDAViewGnuplot2d(da_Stokes_analytic, X_analytic, "Analytic solution for Stokes eqn.", "X_analytic"));
1580 PetscCall(DMDAIntegrateErrors(da_Stokes_analytic, X, X_analytic));
1581
1582 PetscCall(VecAXPY(X_analytic, -1.0, X));
1583 PetscCall(VecGetSize(X_analytic, &N));
1584 N = N / 3;
1585
1586 PetscCall(VecStrideNorm(X_analytic, 0, NORM_1, &nrm1[0]));
1587 PetscCall(VecStrideNorm(X_analytic, 0, NORM_2, &nrm2[0]));
1588 PetscCall(VecStrideNorm(X_analytic, 0, NORM_INFINITY, &nrmI[0]));
1589
1590 PetscCall(VecStrideNorm(X_analytic, 1, NORM_1, &nrm1[1]));
1591 PetscCall(VecStrideNorm(X_analytic, 1, NORM_2, &nrm2[1]));
1592 PetscCall(VecStrideNorm(X_analytic, 1, NORM_INFINITY, &nrmI[1]));
1593
1594 PetscCall(VecStrideNorm(X_analytic, 2, NORM_1, &nrm1[2]));
1595 PetscCall(VecStrideNorm(X_analytic, 2, NORM_2, &nrm2[2]));
1596 PetscCall(VecStrideNorm(X_analytic, 2, NORM_INFINITY, &nrmI[2]));
1597
1598 PetscCall(DMDestroy(&da_Stokes_analytic));
1599 PetscCall(VecDestroy(&X_analytic));
1600 }
1601
1602 PetscCall(KSPDestroy(&ksp_S));
1603 PetscCall(VecDestroy(&X));
1604 PetscCall(VecDestroy(&f));
1605 PetscCall(MatDestroy(&A));
1606 PetscCall(MatDestroy(&B));
1607
1608 PetscCall(DMDestroy(&da_Stokes));
1609 PetscCall(DMDestroy(&da_prop));
1610
1611 PetscCall(VecDestroy(&properties));
1612 PetscCall(VecDestroy(&l_properties));
1613 PetscFunctionReturn(PETSC_SUCCESS);
1614 }
1615
main(int argc,char ** args)1616 int main(int argc, char **args)
1617 {
1618 PetscInt mx, my;
1619
1620 PetscFunctionBeginUser;
1621 PetscCall(PetscInitialize(&argc, &args, NULL, help));
1622 mx = my = 10;
1623 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1624 PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1625 PetscCall(solve_stokes_2d_coupled(mx, my));
1626 PetscCall(PetscFinalize());
1627 return 0;
1628 }
1629
1630 /* -------------------------- helpers for boundary conditions -------------------------------- */
BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)1631 static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
1632 {
1633 DM cda;
1634 Vec coords;
1635 PetscInt si, sj, nx, ny, i, j;
1636 PetscInt M, N;
1637 DMDACoor2d **_coords;
1638 const PetscInt *g_idx;
1639 PetscInt *bc_global_ids;
1640 PetscScalar *bc_vals;
1641 PetscInt nbcs;
1642 PetscInt n_dofs;
1643 ISLocalToGlobalMapping ltogm;
1644
1645 PetscFunctionBeginUser;
1646 PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1647 PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1648
1649 PetscCall(DMGetCoordinateDM(da, &cda));
1650 PetscCall(DMGetCoordinatesLocal(da, &coords));
1651 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1652 PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1653 PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1654
1655 PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1656 PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1657
1658 /* init the entries to -1 so VecSetValues will ignore them */
1659 for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1660
1661 i = nx - 1;
1662 for (j = 0; j < ny; j++) {
1663 PetscInt local_id;
1664
1665 local_id = i + j * nx;
1666
1667 bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1668
1669 bc_vals[j] = 0.0;
1670 }
1671 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1672 nbcs = 0;
1673 if ((si + nx) == (M)) nbcs = ny;
1674
1675 if (b) {
1676 PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1677 PetscCall(VecAssemblyBegin(b));
1678 PetscCall(VecAssemblyEnd(b));
1679 }
1680 if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1681
1682 PetscCall(PetscFree(bc_vals));
1683 PetscCall(PetscFree(bc_global_ids));
1684
1685 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1686 PetscFunctionReturn(PETSC_SUCCESS);
1687 }
1688
BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)1689 static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
1690 {
1691 DM cda;
1692 Vec coords;
1693 PetscInt si, sj, nx, ny, i, j;
1694 PetscInt M, N;
1695 DMDACoor2d **_coords;
1696 const PetscInt *g_idx;
1697 PetscInt *bc_global_ids;
1698 PetscScalar *bc_vals;
1699 PetscInt nbcs;
1700 PetscInt n_dofs;
1701 ISLocalToGlobalMapping ltogm;
1702
1703 PetscFunctionBeginUser;
1704 PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1705 PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1706
1707 PetscCall(DMGetCoordinateDM(da, &cda));
1708 PetscCall(DMGetCoordinatesLocal(da, &coords));
1709 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1710 PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1711 PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1712
1713 PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1714 PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1715
1716 /* init the entries to -1 so VecSetValues will ignore them */
1717 for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1718
1719 i = 0;
1720 for (j = 0; j < ny; j++) {
1721 PetscInt local_id;
1722
1723 local_id = i + j * nx;
1724
1725 bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1726
1727 bc_vals[j] = 0.0;
1728 }
1729 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1730 nbcs = 0;
1731 if (si == 0) nbcs = ny;
1732
1733 if (b) {
1734 PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1735 PetscCall(VecAssemblyBegin(b));
1736 PetscCall(VecAssemblyEnd(b));
1737 }
1738
1739 if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1740
1741 PetscCall(PetscFree(bc_vals));
1742 PetscCall(PetscFree(bc_global_ids));
1743
1744 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1745 PetscFunctionReturn(PETSC_SUCCESS);
1746 }
1747
BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)1748 static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
1749 {
1750 DM cda;
1751 Vec coords;
1752 PetscInt si, sj, nx, ny, i, j;
1753 PetscInt M, N;
1754 DMDACoor2d **_coords;
1755 const PetscInt *g_idx;
1756 PetscInt *bc_global_ids;
1757 PetscScalar *bc_vals;
1758 PetscInt nbcs;
1759 PetscInt n_dofs;
1760 ISLocalToGlobalMapping ltogm;
1761
1762 PetscFunctionBeginUser;
1763 PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1764 PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1765
1766 PetscCall(DMGetCoordinateDM(da, &cda));
1767 PetscCall(DMGetCoordinatesLocal(da, &coords));
1768 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1769 PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1770 PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1771
1772 PetscCall(PetscMalloc1(nx, &bc_global_ids));
1773 PetscCall(PetscMalloc1(nx, &bc_vals));
1774
1775 /* init the entries to -1 so VecSetValues will ignore them */
1776 for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1777
1778 j = ny - 1;
1779 for (i = 0; i < nx; i++) {
1780 PetscInt local_id;
1781
1782 local_id = i + j * nx;
1783
1784 bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1785
1786 bc_vals[i] = 0.0;
1787 }
1788 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1789 nbcs = 0;
1790 if ((sj + ny) == (N)) nbcs = nx;
1791
1792 if (b) {
1793 PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1794 PetscCall(VecAssemblyBegin(b));
1795 PetscCall(VecAssemblyEnd(b));
1796 }
1797 if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1798
1799 PetscCall(PetscFree(bc_vals));
1800 PetscCall(PetscFree(bc_global_ids));
1801
1802 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1803 PetscFunctionReturn(PETSC_SUCCESS);
1804 }
1805
BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)1806 static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
1807 {
1808 DM cda;
1809 Vec coords;
1810 PetscInt si, sj, nx, ny, i, j;
1811 PetscInt M, N;
1812 DMDACoor2d **_coords;
1813 const PetscInt *g_idx;
1814 PetscInt *bc_global_ids;
1815 PetscScalar *bc_vals;
1816 PetscInt nbcs;
1817 PetscInt n_dofs;
1818 ISLocalToGlobalMapping ltogm;
1819
1820 PetscFunctionBeginUser;
1821 PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1822 PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1823
1824 PetscCall(DMGetCoordinateDM(da, &cda));
1825 PetscCall(DMGetCoordinatesLocal(da, &coords));
1826 PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1827 PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1828 PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1829
1830 PetscCall(PetscMalloc1(nx, &bc_global_ids));
1831 PetscCall(PetscMalloc1(nx, &bc_vals));
1832
1833 /* init the entries to -1 so VecSetValues will ignore them */
1834 for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1835
1836 j = 0;
1837 for (i = 0; i < nx; i++) {
1838 PetscInt local_id;
1839
1840 local_id = i + j * nx;
1841
1842 bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1843
1844 bc_vals[i] = 0.0;
1845 }
1846 PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1847 nbcs = 0;
1848 if (sj == 0) nbcs = nx;
1849
1850 if (b) {
1851 PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1852 PetscCall(VecAssemblyBegin(b));
1853 PetscCall(VecAssemblyEnd(b));
1854 }
1855 if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1856
1857 PetscCall(PetscFree(bc_vals));
1858 PetscCall(PetscFree(bc_global_ids));
1859
1860 PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1861 PetscFunctionReturn(PETSC_SUCCESS);
1862 }
1863
1864 /* Impose free slip boundary conditions; u_{i} n_{i} = 0, tau_{ij} t_j = 0 */
DMDABCApplyFreeSlip(DM da_Stokes,Mat A,Vec f)1865 static PetscErrorCode DMDABCApplyFreeSlip(DM da_Stokes, Mat A, Vec f)
1866 {
1867 PetscFunctionBeginUser;
1868 PetscCall(BCApplyZero_NORTH(da_Stokes, 1, A, f));
1869 PetscCall(BCApplyZero_EAST(da_Stokes, 0, A, f));
1870 PetscCall(BCApplyZero_SOUTH(da_Stokes, 1, A, f));
1871 PetscCall(BCApplyZero_WEST(da_Stokes, 0, A, f));
1872 PetscFunctionReturn(PETSC_SUCCESS);
1873 }
1874
1875 /*TEST
1876
1877 build:
1878 requires: !complex !single
1879
1880 test:
1881 args: -stokes_pc_use_amat -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type lu -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short
1882
1883 testset:
1884 args: -stokes_pc_use_amat -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_fieldsplit_u_ksp_type preonly -stokes_fieldsplit_u_pc_type lu -stokes_fieldsplit_p_ksp_type preonly -stokes_fieldsplit_p_pc_type jacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short
1885 test:
1886 suffix: 2
1887 args:
1888 output_file: output/ex43_1.out
1889 test:
1890 requires: mumps
1891 suffix: 2_mumps
1892 args: -change true -stokes_ksp_view
1893 output_file: output/ex43_2_mumps.out
1894 # mumps INFO,INFOG,RINFO,RINFOG may vary on different archs, so keep just a stable one
1895 filter: grep -E -v "(INFOG\([^7]|INFO\(|\[0\])" | grep -v " total: nonzeros="
1896
1897 test:
1898 suffix: 3
1899 nsize: 4
1900 args: -stokes_pc_use_amat -stokes_ksp_type gcr -stokes_ksp_gcr_restart 60 -stokes_ksp_norm_type unpreconditioned -stokes_ksp_rtol 1.e-2 -c_str 3 -sinker_eta0 1.0 -sinker_eta1 100 -sinker_dx 0.4 -sinker_dy 0.3 -mx 128 -my 128 -stokes_ksp_monitor_short -stokes_pc_type mg -stokes_mg_levels_pc_type fieldsplit -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_levels_pc_fieldsplit_block_size 3 -stokes_mg_levels_pc_fieldsplit_0_fields 0,1 -stokes_mg_levels_pc_fieldsplit_1_fields 2 -stokes_mg_levels_fieldsplit_0_pc_type sor -stokes_mg_levels_fieldsplit_1_pc_type sor -stokes_mg_levels_ksp_type chebyshev -stokes_mg_levels_ksp_max_it 1 -stokes_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -stokes_pc_mg_levels 4 -stokes_ksp_view
1901
1902 test:
1903 suffix: 4
1904 nsize: 4
1905 args: -stokes_ksp_type pipegcr -stokes_ksp_pipegcr_mmax 60 -stokes_ksp_pipegcr_unroll_w 1 -stokes_ksp_norm_type natural -c_str 3 -sinker_eta0 1.0 -sinker_eta1 100 -sinker_dx 0.4 -sinker_dy 0.3 -mx 128 -my 128 -stokes_ksp_monitor_short -stokes_pc_type mg -stokes_mg_levels_pc_type fieldsplit -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_levels_pc_fieldsplit_block_size 3 -stokes_mg_levels_pc_fieldsplit_0_fields 0,1 -stokes_mg_levels_pc_fieldsplit_1_fields 2 -stokes_mg_levels_fieldsplit_0_pc_type sor -stokes_mg_levels_fieldsplit_1_pc_type sor -stokes_mg_levels_ksp_type chebyshev -stokes_mg_levels_ksp_max_it 1 -stokes_mg_levels_ksp_chebyshev_esteig 0,0.2,0,1.1 -stokes_pc_mg_levels 4 -stokes_ksp_view
1906
1907 test:
1908 suffix: 5
1909 nsize: 4
1910 args: -stokes_pc_fieldsplit_off_diag_use_amat -stokes_ksp_type pipegcr -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type bjacobi -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type bjacobi -c_str 0 -solcx_eta0 1.0 -solcx_eta1 1.0e6 -solcx_xc 0.5 -solcx_nz 2 -mx 20 -my 20 -stokes_ksp_monitor_short -stokes_ksp_view
1911
1912 test:
1913 suffix: 6
1914 nsize: 8
1915 args: -stokes_ksp_view -stokes_pc_type mg -stokes_pc_mg_levels 2 -stokes_mg_coarse_pc_type telescope -stokes_mg_coarse_pc_telescope_reduction_factor 2 -stokes_pc_use_amat false -stokes_pc_mg_galerkin pmat -stokes_mg_coarse_pc_telescope_subcomm_type contiguous
1916
1917 test:
1918 suffix: bjacobi
1919 nsize: 4
1920 args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type aij -stokes_ksp_converged_reason
1921
1922 test:
1923 suffix: bjacobi_baij
1924 nsize: 4
1925 args: -stokes_ksp_rtol 1.e-4 -stokes_pc_type bjacobi -stokes_pc_bjacobi_blocks 2 -dm_mat_type baij -stokes_ksp_converged_reason
1926 output_file: output/ex43_bjacobi.out
1927
1928 test:
1929 suffix: nested_gmg
1930 nsize: 4
1931 args: -stokes_pc_fieldsplit_off_diag_use_amat -mx 16 -my 16 -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_fieldsplit_u_pc_type mg -stokes_fieldsplit_u_pc_mg_levels 5 -stokes_fieldsplit_u_pc_mg_galerkin pmat -stokes_fieldsplit_u_ksp_type cg -stokes_fieldsplit_u_ksp_rtol 1.0e-4 -stokes_fieldsplit_u_mg_levels_pc_type jacobi -solcx_eta0 1.0e4 -stokes_fieldsplit_u_ksp_converged_reason -stokes_ksp_converged_reason -stokes_fieldsplit_p_sub_pc_factor_zeropivot 1.e-8
1932
1933 test:
1934 suffix: fetidp
1935 nsize: 8
1936 args: -dm_mat_type is -stokes_ksp_type fetidp -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_converged_reason -stokes_fetidp_pc_fieldsplit_schur_fact_type diag -stokes_fetidp_fieldsplit_p_pc_type bjacobi -stokes_fetidp_fieldsplit_lag_ksp_type preonly -stokes_fetidp_fieldsplit_p_ksp_type preonly -stokes_ksp_fetidp_pressure_field 2 -stokes_fetidp_pc_fieldsplit_schur_scale -1
1937
1938 test:
1939 suffix: fetidp_unsym
1940 nsize: 8
1941 args: -dm_mat_type is -stokes_ksp_type fetidp -stokes_ksp_monitor_true_residual -stokes_ksp_converged_reason -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
1942
1943 test:
1944 suffix: bddc_stokes_deluxe
1945 nsize: 8
1946 args: -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -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
1947
1948 test:
1949 suffix: bddc_stokes_subdomainjump_deluxe
1950 nsize: 9
1951 args: -c_str 4 -jump_magnitude 3 -stokes_ksp_monitor_short -stokes_ksp_converged_reason -stokes_pc_type bddc -dm_mat_type is -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
1952
1953 TEST*/
1954