Actual source code: ex70.c
1: static char help[] = "------------------------------------------------------------------------------------------------------------------------------ \n\
2: Solves the time-dependent incompressible, variable viscosity Stokes equation in 2D driven by buouyancy variations. \n\
3: Time-dependence is introduced by evolving the density (rho) and viscosity (eta) according to \n\
4: D \\rho / Dt = 0 and D \\eta / Dt = 0 \n\
5: The Stokes problem is discretized using Q1-Q1 finite elements, stabilized with Bochev's polynomial projection method. \n\
6: The hyperbolic evolution equation for density is discretized using a variant of the Particle-In-Cell (PIC) method. \n\
7: The DMDA object is used to define the FE problem, whilst DMSwarm provides support for the PIC method. \n\
8: Material points (particles) store density and viscosity. The particles are advected with the fluid velocity using RK1. \n\
9: At each time step, the value of density and viscosity stored on each particle are projected into a Q1 function space \n\
10: and then interpolated onto the Gauss quadrature points. \n\
11: The model problem defined in this example is the iso-viscous Rayleigh-Taylor instability (case 1a) from: \n\
12: \"A comparison of methods for the modeling of thermochemical convection\" \n\
13: P.E. van Keken, S.D. King, H. Schmeling, U.R. Christensen, D. Neumeister and M.-P. Doin, \n\
14: Journal of Geophysical Research, vol 102 (B10), 477--499 (1997) \n\
15: Note that whilst the model problem defined is for an iso-viscous, the implementation in this example supports \n\
16: variable viscoity formulations. \n\
17: This example is based on src/ksp/ksp/tutorials/ex43.c \n\
18: Options: \n\
19: -mx : Number of elements in the x-direction \n\
20: -my : Number of elements in the y-direction \n\
21: -mxy : Number of elements in the x- and y-directions \n\
22: -nt : Number of time steps \n\
23: -dump_freq : Frequency of output file creation \n\
24: -ppcell : Number of times the reference cell is sub-divided \n\
25: -randomize_coords : Apply a random shift to each particle coordinate in the range [-fac*dh,0.fac*dh] \n\
26: -randomize_fac : Set the scaling factor for the random shift (default = 0.25)\n";
28: /* Contributed by Dave May */
30: #include <petscksp.h>
31: #include <petscdm.h>
32: #include <petscdmda.h>
33: #include <petscdmswarm.h>
34: #include <petsc/private/dmimpl.h>
36: static PetscErrorCode DMDAApplyBoundaryConditions(DM, Mat, Vec);
38: #define NSD 2 /* number of spatial dimensions */
39: #define NODES_PER_EL 4 /* nodes per element */
40: #define U_DOFS 2 /* degrees of freedom per velocity node */
41: #define P_DOFS 1 /* degrees of freedom per pressure node */
42: #define GAUSS_POINTS 4
44: static void EvaluateBasis_Q1(PetscScalar _xi[], PetscScalar N[])
45: {
46: PetscScalar xi = _xi[0];
47: PetscScalar eta = _xi[1];
49: N[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
50: N[1] = 0.25 * (1.0 + xi) * (1.0 - eta);
51: N[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
52: N[3] = 0.25 * (1.0 - xi) * (1.0 + eta);
53: }
55: static void EvaluateBasisDerivatives_Q1(PetscScalar _xi[], PetscScalar dN[][NODES_PER_EL])
56: {
57: PetscScalar xi = _xi[0];
58: PetscScalar eta = _xi[1];
60: dN[0][0] = -0.25 * (1.0 - eta);
61: dN[0][1] = 0.25 * (1.0 - eta);
62: dN[0][2] = 0.25 * (1.0 + eta);
63: dN[0][3] = -0.25 * (1.0 + eta);
65: dN[1][0] = -0.25 * (1.0 - xi);
66: dN[1][1] = -0.25 * (1.0 + xi);
67: dN[1][2] = 0.25 * (1.0 + xi);
68: dN[1][3] = 0.25 * (1.0 - xi);
69: }
71: static void EvaluateDerivatives(PetscScalar dN[][NODES_PER_EL], PetscScalar dNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
72: {
73: PetscScalar J00, J01, J10, J11, J;
74: PetscScalar iJ00, iJ01, iJ10, iJ11;
75: PetscInt i;
77: J00 = J01 = J10 = J11 = 0.0;
78: for (i = 0; i < NODES_PER_EL; i++) {
79: PetscScalar cx = coords[2 * i];
80: PetscScalar cy = coords[2 * i + 1];
82: J00 += dN[0][i] * cx; /* J_xx = dx/dxi */
83: J01 += dN[0][i] * cy; /* J_xy = dy/dxi */
84: J10 += dN[1][i] * cx; /* J_yx = dx/deta */
85: J11 += dN[1][i] * cy; /* J_yy = dy/deta */
86: }
87: J = (J00 * J11) - (J01 * J10);
89: iJ00 = J11 / J;
90: iJ01 = -J01 / J;
91: iJ10 = -J10 / J;
92: iJ11 = J00 / J;
94: for (i = 0; i < NODES_PER_EL; i++) {
95: dNx[0][i] = dN[0][i] * iJ00 + dN[1][i] * iJ01;
96: dNx[1][i] = dN[0][i] * iJ10 + dN[1][i] * iJ11;
97: }
99: if (det_J) *det_J = J;
100: }
102: static void CreateGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
103: {
104: *ngp = 4;
105: gp_xi[0][0] = -0.57735026919;
106: gp_xi[0][1] = -0.57735026919;
107: gp_xi[1][0] = -0.57735026919;
108: gp_xi[1][1] = 0.57735026919;
109: gp_xi[2][0] = 0.57735026919;
110: gp_xi[2][1] = 0.57735026919;
111: gp_xi[3][0] = 0.57735026919;
112: gp_xi[3][1] = -0.57735026919;
113: gp_weight[0] = 1.0;
114: gp_weight[1] = 1.0;
115: gp_weight[2] = 1.0;
116: gp_weight[3] = 1.0;
117: }
119: static PetscErrorCode DMDAGetElementEqnums_up(const PetscInt element[], PetscInt s_u[], PetscInt s_p[])
120: {
121: PetscInt i;
122: PetscFunctionBeginUser;
123: for (i = 0; i < NODES_PER_EL; i++) {
124: /* velocity */
125: s_u[NSD * i + 0] = 3 * element[i];
126: s_u[NSD * i + 1] = 3 * element[i] + 1;
127: /* pressure */
128: s_p[i] = 3 * element[i] + 2;
129: }
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: static PetscInt map_wIwDI_uJuDJ(PetscInt wi, PetscInt wd, PetscInt w_NPE, PetscInt w_dof, PetscInt ui, PetscInt ud, PetscInt u_NPE, PetscInt u_dof)
134: {
135: PetscInt ij, r, c, nc;
137: nc = u_NPE * u_dof;
138: r = w_dof * wi + wd;
139: c = u_dof * ui + ud;
140: ij = r * nc + c;
141: return (ij);
142: }
144: static void BForm_DivT(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
145: {
146: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
147: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
148: PetscScalar J_p, tildeD[3];
149: PetscScalar B[3][U_DOFS * NODES_PER_EL];
150: PetscInt p, i, j, k, ngp;
152: /* define quadrature rule */
153: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
155: /* evaluate bilinear form */
156: for (p = 0; p < ngp; p++) {
157: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
158: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
160: for (i = 0; i < NODES_PER_EL; i++) {
161: PetscScalar d_dx_i = GNx_p[0][i];
162: PetscScalar d_dy_i = GNx_p[1][i];
164: B[0][2 * i] = d_dx_i;
165: B[0][2 * i + 1] = 0.0;
166: B[1][2 * i] = 0.0;
167: B[1][2 * i + 1] = d_dy_i;
168: B[2][2 * i] = d_dy_i;
169: B[2][2 * i + 1] = d_dx_i;
170: }
172: tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
173: tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
174: tildeD[2] = gp_weight[p] * J_p * eta[p];
176: /* form Bt tildeD B */
177: /*
178: Ke_ij = Bt_ik . D_kl . B_lj
179: = B_ki . D_kl . B_lj
180: = B_ki . D_kk . B_kj
181: */
182: for (i = 0; i < 8; i++) {
183: for (j = 0; j < 8; j++) {
184: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
185: Ke[i + 8 * j] += B[k][i] * tildeD[k] * B[k][j];
186: }
187: }
188: }
189: }
190: }
192: static void BForm_Grad(PetscScalar Ke[], PetscScalar coords[])
193: {
194: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
195: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
196: PetscScalar J_p, fac;
197: PetscInt p, i, j, di, ngp;
199: /* define quadrature rule */
200: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
202: /* evaluate bilinear form */
203: for (p = 0; p < ngp; p++) {
204: EvaluateBasis_Q1(gp_xi[p], Ni_p);
205: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
206: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
207: fac = gp_weight[p] * J_p;
209: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
210: for (di = 0; di < NSD; di++) { /* u dofs */
211: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
212: PetscInt IJ;
213: IJ = map_wIwDI_uJuDJ(i, di, NODES_PER_EL, 2, j, 0, NODES_PER_EL, 1);
215: Ke[IJ] -= GNx_p[di][i] * Ni_p[j] * fac;
216: }
217: }
218: }
219: }
220: }
222: static void BForm_Div(PetscScalar De[], PetscScalar coords[])
223: {
224: PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
225: PetscInt i, j, nr_g, nc_g;
227: PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));
228: BForm_Grad(Ge, coords);
230: nr_g = U_DOFS * NODES_PER_EL;
231: nc_g = P_DOFS * NODES_PER_EL;
233: for (i = 0; i < nr_g; i++) {
234: for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
235: }
236: }
238: static void BForm_Stabilisation(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
239: {
240: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
241: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
242: PetscScalar J_p, fac, eta_avg;
243: PetscInt p, i, j, ngp;
245: /* define quadrature rule */
246: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
248: /* evaluate bilinear form */
249: for (p = 0; p < ngp; p++) {
250: EvaluateBasis_Q1(gp_xi[p], Ni_p);
251: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
252: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
253: fac = gp_weight[p] * J_p;
255: for (i = 0; i < NODES_PER_EL; i++) {
256: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * (Ni_p[i] * Ni_p[j] - 0.0625);
257: }
258: }
260: /* scale */
261: eta_avg = 0.0;
262: for (p = 0; p < ngp; p++) eta_avg += eta[p];
263: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
264: fac = 1.0 / eta_avg;
265: for (i = 0; i < NODES_PER_EL; i++) {
266: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
267: }
268: }
270: static void BForm_ScaledMassMatrix(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
271: {
272: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
273: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
274: PetscScalar J_p, fac, eta_avg;
275: PetscInt p, i, j, ngp;
277: /* define quadrature rule */
278: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
280: /* evaluate bilinear form */
281: for (p = 0; p < ngp; p++) {
282: EvaluateBasis_Q1(gp_xi[p], Ni_p);
283: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
284: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
285: fac = gp_weight[p] * J_p;
287: for (i = 0; i < NODES_PER_EL; i++) {
288: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * Ni_p[i] * Ni_p[j];
289: }
290: }
292: /* scale */
293: eta_avg = 0.0;
294: for (p = 0; p < ngp; p++) eta_avg += eta[p];
295: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
296: fac = 1.0 / eta_avg;
297: for (i = 0; i < NODES_PER_EL; i++) {
298: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] *= fac;
299: }
300: }
302: static void LForm_MomentumRHS(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
303: {
304: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
305: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
306: PetscScalar J_p, fac;
307: PetscInt p, i, ngp;
309: /* define quadrature rule */
310: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
312: /* evaluate linear form */
313: for (p = 0; p < ngp; p++) {
314: EvaluateBasis_Q1(gp_xi[p], Ni_p);
315: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
316: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
317: fac = gp_weight[p] * J_p;
319: for (i = 0; i < NODES_PER_EL; i++) {
320: Fe[NSD * i] = 0.0;
321: Fe[NSD * i + 1] -= fac * Ni_p[i] * fy[p];
322: }
323: }
324: }
326: static PetscErrorCode GetElementCoords(const PetscScalar _coords[], const PetscInt e2n[], PetscScalar el_coords[])
327: {
328: PetscInt i, d;
329: PetscFunctionBeginUser;
330: /* get coords for the element */
331: for (i = 0; i < 4; i++) {
332: for (d = 0; d < NSD; d++) el_coords[NSD * i + d] = _coords[NSD * e2n[i] + d];
333: }
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: static PetscErrorCode AssembleStokes_A(Mat A, DM stokes_da, DM quadrature)
338: {
339: DM cda;
340: Vec coords;
341: const PetscScalar *_coords;
342: PetscInt u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
343: PetscInt p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
344: PetscInt nel, npe, eidx;
345: const PetscInt *element_list;
346: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
347: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
348: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
349: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
350: PetscScalar el_coords[NODES_PER_EL * NSD];
351: PetscScalar *q_eta, *prop_eta;
353: PetscFunctionBeginUser;
354: PetscCall(MatZeroEntries(A));
355: /* setup for coords */
356: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
357: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
358: PetscCall(VecGetArrayRead(coords, &_coords));
360: /* setup for coefficients */
361: PetscCall(DMSwarmGetField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
363: PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
364: for (eidx = 0; eidx < nel; eidx++) {
365: const PetscInt *element = &element_list[npe * eidx];
367: /* get coords for the element */
368: PetscCall(GetElementCoords(_coords, element, el_coords));
370: /* get coefficients for the element */
371: prop_eta = &q_eta[GAUSS_POINTS * eidx];
373: /* initialise element stiffness matrix */
374: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
375: PetscCall(PetscMemzero(Ge, sizeof(Ge)));
376: PetscCall(PetscMemzero(De, sizeof(De)));
377: PetscCall(PetscMemzero(Ce, sizeof(Ce)));
379: /* form element stiffness matrix */
380: BForm_DivT(Ae, el_coords, prop_eta);
381: BForm_Grad(Ge, el_coords);
382: BForm_Div(De, el_coords);
383: BForm_Stabilisation(Ce, el_coords, prop_eta);
385: /* insert element matrix into global matrix */
386: PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
387: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
388: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
389: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
390: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
391: }
392: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
393: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
395: PetscCall(DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
396: PetscCall(VecRestoreArrayRead(coords, &_coords));
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: static PetscErrorCode AssembleStokes_PC(Mat A, DM stokes_da, DM quadrature)
401: {
402: DM cda;
403: Vec coords;
404: const PetscScalar *_coords;
405: PetscInt u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
406: PetscInt p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
407: PetscInt nel, npe, eidx;
408: const PetscInt *element_list;
409: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
410: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
411: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
412: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
413: PetscScalar el_coords[NODES_PER_EL * NSD];
414: PetscScalar *q_eta, *prop_eta;
416: PetscFunctionBeginUser;
417: PetscCall(MatZeroEntries(A));
418: /* setup for coords */
419: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
420: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
421: PetscCall(VecGetArrayRead(coords, &_coords));
423: /* setup for coefficients */
424: PetscCall(DMSwarmGetField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
426: PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
427: for (eidx = 0; eidx < nel; eidx++) {
428: const PetscInt *element = &element_list[npe * eidx];
430: /* get coords for the element */
431: PetscCall(GetElementCoords(_coords, element, el_coords));
433: /* get coefficients for the element */
434: prop_eta = &q_eta[GAUSS_POINTS * eidx];
436: /* initialise element stiffness matrix */
437: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
438: PetscCall(PetscMemzero(Ge, sizeof(Ge)));
439: PetscCall(PetscMemzero(De, sizeof(De)));
440: PetscCall(PetscMemzero(Ce, sizeof(Ce)));
442: /* form element stiffness matrix */
443: BForm_DivT(Ae, el_coords, prop_eta);
444: BForm_Grad(Ge, el_coords);
445: BForm_ScaledMassMatrix(Ce, el_coords, prop_eta);
447: /* insert element matrix into global matrix */
448: PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
449: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
450: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
451: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
452: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
453: }
454: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
455: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
457: PetscCall(DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
458: PetscCall(VecRestoreArrayRead(coords, &_coords));
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: static PetscErrorCode AssembleStokes_RHS(Vec F, DM stokes_da, DM quadrature)
464: {
465: DM cda;
466: Vec coords;
467: const PetscScalar *_coords;
468: PetscInt u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
469: PetscInt p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
470: PetscInt nel, npe, eidx, i;
471: const PetscInt *element_list;
472: PetscScalar Fe[NODES_PER_EL * U_DOFS];
473: PetscScalar He[NODES_PER_EL * P_DOFS];
474: PetscScalar el_coords[NODES_PER_EL * NSD];
475: PetscScalar *q_rhs, *prop_fy;
476: Vec local_F;
477: PetscScalar *LA_F;
479: PetscFunctionBeginUser;
480: PetscCall(VecZeroEntries(F));
481: /* setup for coords */
482: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
483: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
484: PetscCall(VecGetArrayRead(coords, &_coords));
486: /* setup for coefficients */
487: PetscCall(DMSwarmGetField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
489: /* get access to the vector */
490: PetscCall(DMGetLocalVector(stokes_da, &local_F));
491: PetscCall(VecZeroEntries(local_F));
492: PetscCall(VecGetArray(local_F, &LA_F));
494: PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
495: for (eidx = 0; eidx < nel; eidx++) {
496: const PetscInt *element = &element_list[npe * eidx];
498: /* get coords for the element */
499: PetscCall(GetElementCoords(_coords, element, el_coords));
501: /* get coefficients for the element */
502: prop_fy = &q_rhs[GAUSS_POINTS * eidx];
504: /* initialise element stiffness matrix */
505: PetscCall(PetscMemzero(Fe, sizeof(Fe)));
506: PetscCall(PetscMemzero(He, sizeof(He)));
508: /* form element stiffness matrix */
509: LForm_MomentumRHS(Fe, el_coords, NULL, prop_fy);
511: /* insert element matrix into global matrix */
512: PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
514: for (i = 0; i < NODES_PER_EL * U_DOFS; i++) LA_F[u_eqn[i]] += Fe[i];
515: }
516: PetscCall(DMSwarmRestoreField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
517: PetscCall(VecRestoreArrayRead(coords, &_coords));
519: PetscCall(VecRestoreArray(local_F, &LA_F));
520: PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
521: PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
522: PetscCall(DMRestoreLocalVector(stokes_da, &local_F));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: PetscErrorCode DMSwarmPICInsertPointsCellwise(DM dm, DM dmc, PetscInt e, PetscInt npoints, PetscReal xi[], PetscBool proximity_initialization)
528: {
529: PetscInt dim, nel, npe, q, k, d, ncurr;
530: const PetscInt *element_list;
531: Vec coor;
532: const PetscScalar *_coor;
533: PetscReal **basis, *elcoor, *xp;
534: PetscReal *swarm_coor;
535: PetscInt *swarm_cellid;
537: PetscFunctionBeginUser;
538: PetscCall(DMGetDimension(dm, &dim));
539: PetscCall(DMDAGetElements(dmc, &nel, &npe, &element_list));
541: PetscCall(PetscMalloc1(dim * npoints, &xp));
542: PetscCall(PetscMalloc1(dim * npe, &elcoor));
543: PetscCall(PetscMalloc1(npoints, &basis));
544: for (q = 0; q < npoints; q++) {
545: PetscCall(PetscMalloc1(npe, &basis[q]));
547: switch (dim) {
548: case 1:
549: basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]);
550: basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]);
551: break;
552: case 2:
553: basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
554: basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
555: basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
556: basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
557: break;
559: case 3:
560: basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
561: basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
562: basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
563: basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
564: basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
565: basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
566: basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
567: basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
568: break;
569: }
570: }
572: PetscCall(DMGetCoordinatesLocal(dmc, &coor));
573: PetscCall(VecGetArrayRead(coor, &_coor));
574: /* compute and store the coordinates for the new points */
575: {
576: const PetscInt *element = &element_list[npe * e];
578: for (k = 0; k < npe; k++) {
579: for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]);
580: }
581: for (q = 0; q < npoints; q++) {
582: for (d = 0; d < dim; d++) xp[dim * q + d] = 0.0;
583: for (k = 0; k < npe; k++) {
584: for (d = 0; d < dim; d++) xp[dim * q + d] += basis[q][k] * elcoor[dim * k + d];
585: }
586: }
587: }
588: PetscCall(VecRestoreArrayRead(coor, &_coor));
589: PetscCall(DMDARestoreElements(dmc, &nel, &npe, &element_list));
591: PetscCall(DMSwarmGetLocalSize(dm, &ncurr));
592: PetscCall(DMSwarmAddNPoints(dm, npoints));
594: if (proximity_initialization) {
595: PetscInt *nnlist;
596: PetscReal *coor_q, *coor_qn;
597: PetscInt npoints_e, *plist_e;
599: PetscCall(DMSwarmSortGetPointsPerCell(dm, e, &npoints_e, &plist_e));
601: PetscCall(PetscMalloc1(npoints, &nnlist));
602: /* find nearest neighbour points in this cell */
603: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
604: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
605: for (q = 0; q < npoints; q++) {
606: PetscInt qn, nearest_neighbour = -1;
607: PetscReal sep, min_sep = PETSC_MAX_REAL;
609: coor_q = &xp[dim * q];
610: for (qn = 0; qn < npoints_e; qn++) {
611: coor_qn = &swarm_coor[dim * plist_e[qn]];
612: sep = 0.0;
613: for (d = 0; d < dim; d++) sep += (coor_q[d] - coor_qn[d]) * (coor_q[d] - coor_qn[d]);
614: if (sep < min_sep) {
615: nearest_neighbour = plist_e[qn];
616: min_sep = sep;
617: }
618: }
619: PetscCheck(nearest_neighbour != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell %" PetscInt_FMT " is empty - cannot initialize using nearest neighbours", e);
620: nnlist[q] = nearest_neighbour;
621: }
622: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
623: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
625: /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
626: for (q = 0; q < npoints; q++) PetscCall(DMSwarmCopyPoint(dm, nnlist[q], ncurr + q));
627: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
628: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
629: for (q = 0; q < npoints; q++) {
630: /* set the coordinates */
631: for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
632: /* set the cell index */
633: swarm_cellid[ncurr + q] = e;
634: }
635: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
636: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
638: PetscCall(PetscFree(plist_e));
639: PetscCall(PetscFree(nnlist));
640: } else {
641: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
642: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
643: for (q = 0; q < npoints; q++) {
644: /* set the coordinates */
645: for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
646: /* set the cell index */
647: swarm_cellid[ncurr + q] = e;
648: }
649: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
650: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
651: }
653: PetscCall(PetscFree(xp));
654: PetscCall(PetscFree(elcoor));
655: for (q = 0; q < npoints; q++) PetscCall(PetscFree(basis[q]));
656: PetscCall(PetscFree(basis));
657: PetscFunctionReturn(PETSC_SUCCESS);
658: }
660: PetscErrorCode MaterialPoint_PopulateCell(DM dm_vp, DM dm_mpoint)
661: {
662: PetscInt _npe, _nel, e, nel;
663: const PetscInt *element;
664: DM dmc;
665: PetscQuadrature quadrature;
666: const PetscReal *xi;
667: PetscInt npoints_q, cnt, cnt_g;
669: PetscFunctionBeginUser;
670: PetscCall(DMDAGetElements(dm_vp, &_nel, &_npe, &element));
671: nel = _nel;
672: PetscCall(DMDARestoreElements(dm_vp, &_nel, &_npe, &element));
674: PetscCall(PetscDTGaussTensorQuadrature(2, 1, 4, -1.0, 1.0, &quadrature));
675: PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xi, NULL));
676: PetscCall(DMSwarmGetCellDM(dm_mpoint, &dmc));
678: PetscCall(DMSwarmSortGetAccess(dm_mpoint));
680: cnt = 0;
681: for (e = 0; e < nel; e++) {
682: PetscInt npoints_per_cell;
684: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm_mpoint, e, &npoints_per_cell));
686: if (npoints_per_cell < 12) {
687: PetscCall(DMSwarmPICInsertPointsCellwise(dm_mpoint, dm_vp, e, npoints_q, (PetscReal *)xi, PETSC_TRUE));
688: cnt++;
689: }
690: }
691: PetscCallMPI(MPI_Allreduce(&cnt, &cnt_g, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
692: if (cnt_g > 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... ....pop cont: adjusted %" PetscInt_FMT " cells\n", cnt_g));
694: PetscCall(DMSwarmSortRestoreAccess(dm_mpoint));
695: PetscCall(PetscQuadratureDestroy(&quadrature));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: PetscErrorCode MaterialPoint_AdvectRK1(DM dm_vp, Vec vp, PetscReal dt, DM dm_mpoint)
700: {
701: Vec vp_l, coor_l;
702: const PetscScalar *LA_vp;
703: PetscInt i, p, e, npoints, nel, npe;
704: PetscInt *mpfield_cell;
705: PetscReal *mpfield_coor;
706: const PetscInt *element_list;
707: const PetscInt *element;
708: PetscScalar xi_p[NSD], Ni[NODES_PER_EL];
709: const PetscScalar *LA_coor;
710: PetscScalar dx[NSD];
712: PetscFunctionBeginUser;
713: PetscCall(DMGetCoordinatesLocal(dm_vp, &coor_l));
714: PetscCall(VecGetArrayRead(coor_l, &LA_coor));
716: PetscCall(DMGetLocalVector(dm_vp, &vp_l));
717: PetscCall(DMGlobalToLocalBegin(dm_vp, vp, INSERT_VALUES, vp_l));
718: PetscCall(DMGlobalToLocalEnd(dm_vp, vp, INSERT_VALUES, vp_l));
719: PetscCall(VecGetArrayRead(vp_l, &LA_vp));
721: PetscCall(DMDAGetElements(dm_vp, &nel, &npe, &element_list));
722: PetscCall(DMSwarmGetLocalSize(dm_mpoint, &npoints));
723: PetscCall(DMSwarmGetField(dm_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
724: PetscCall(DMSwarmGetField(dm_mpoint, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
725: for (p = 0; p < npoints; p++) {
726: PetscReal *coor_p;
727: PetscScalar vel_n[NSD * NODES_PER_EL], vel_p[NSD];
728: const PetscScalar *x0;
729: const PetscScalar *x2;
731: e = mpfield_cell[p];
732: coor_p = &mpfield_coor[NSD * p];
733: element = &element_list[NODES_PER_EL * e];
735: /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
736: x0 = &LA_coor[NSD * element[0]];
737: x2 = &LA_coor[NSD * element[2]];
739: dx[0] = x2[0] - x0[0];
740: dx[1] = x2[1] - x0[1];
742: xi_p[0] = 2.0 * (coor_p[0] - x0[0]) / dx[0] - 1.0;
743: xi_p[1] = 2.0 * (coor_p[1] - x0[1]) / dx[1] - 1.0;
744: PetscCheck(PetscRealPart(xi_p[0]) >= -1.0 - PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (xi) too small %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[0]), e);
745: PetscCheck(PetscRealPart(xi_p[0]) <= 1.0 + PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (xi) too large %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[0]), e);
746: PetscCheck(PetscRealPart(xi_p[1]) >= -1.0 - PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (eta) too small %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[1]), e);
747: PetscCheck(PetscRealPart(xi_p[1]) <= 1.0 + PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (eta) too large %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[1]), e);
749: /* evaluate basis functions */
750: EvaluateBasis_Q1(xi_p, Ni);
752: /* get cell nodal velocities */
753: for (i = 0; i < NODES_PER_EL; i++) {
754: PetscInt nid;
756: nid = element[i];
757: vel_n[NSD * i + 0] = LA_vp[(NSD + 1) * nid + 0];
758: vel_n[NSD * i + 1] = LA_vp[(NSD + 1) * nid + 1];
759: }
761: /* interpolate velocity */
762: vel_p[0] = vel_p[1] = 0.0;
763: for (i = 0; i < NODES_PER_EL; i++) {
764: vel_p[0] += Ni[i] * vel_n[NSD * i + 0];
765: vel_p[1] += Ni[i] * vel_n[NSD * i + 1];
766: }
768: coor_p[0] += dt * PetscRealPart(vel_p[0]);
769: coor_p[1] += dt * PetscRealPart(vel_p[1]);
770: }
772: PetscCall(DMSwarmRestoreField(dm_mpoint, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
773: PetscCall(DMSwarmRestoreField(dm_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
774: PetscCall(DMDARestoreElements(dm_vp, &nel, &npe, &element_list));
775: PetscCall(VecRestoreArrayRead(vp_l, &LA_vp));
776: PetscCall(DMRestoreLocalVector(dm_vp, &vp_l));
777: PetscCall(VecRestoreArrayRead(coor_l, &LA_coor));
778: PetscFunctionReturn(PETSC_SUCCESS);
779: }
781: PetscErrorCode MaterialPoint_Interpolate(DM dm, Vec eta_v, Vec rho_v, DM dm_quadrature)
782: {
783: Vec eta_l, rho_l;
784: PetscScalar *_eta_l, *_rho_l;
785: PetscInt nqp, npe, nel;
786: PetscScalar qp_xi[GAUSS_POINTS][NSD];
787: PetscScalar qp_weight[GAUSS_POINTS];
788: PetscInt q, k, e;
789: PetscScalar Ni[GAUSS_POINTS][NODES_PER_EL];
790: const PetscInt *element_list;
791: PetscReal *q_eta, *q_rhs;
793: PetscFunctionBeginUser;
794: /* define quadrature rule */
795: CreateGaussQuadrature(&nqp, qp_xi, qp_weight);
796: for (q = 0; q < nqp; q++) EvaluateBasis_Q1(qp_xi[q], Ni[q]);
798: PetscCall(DMGetLocalVector(dm, &eta_l));
799: PetscCall(DMGetLocalVector(dm, &rho_l));
801: PetscCall(DMGlobalToLocalBegin(dm, eta_v, INSERT_VALUES, eta_l));
802: PetscCall(DMGlobalToLocalEnd(dm, eta_v, INSERT_VALUES, eta_l));
803: PetscCall(DMGlobalToLocalBegin(dm, rho_v, INSERT_VALUES, rho_l));
804: PetscCall(DMGlobalToLocalEnd(dm, rho_v, INSERT_VALUES, rho_l));
806: PetscCall(VecGetArray(eta_l, &_eta_l));
807: PetscCall(VecGetArray(rho_l, &_rho_l));
809: PetscCall(DMSwarmGetField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
810: PetscCall(DMSwarmGetField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
812: PetscCall(DMDAGetElements(dm, &nel, &npe, &element_list));
813: for (e = 0; e < nel; e++) {
814: PetscScalar eta_field_e[NODES_PER_EL];
815: PetscScalar rho_field_e[NODES_PER_EL];
816: const PetscInt *element = &element_list[4 * e];
818: for (k = 0; k < NODES_PER_EL; k++) {
819: eta_field_e[k] = _eta_l[element[k]];
820: rho_field_e[k] = _rho_l[element[k]];
821: }
823: for (q = 0; q < nqp; q++) {
824: PetscScalar eta_q, rho_q;
826: eta_q = rho_q = 0.0;
827: for (k = 0; k < NODES_PER_EL; k++) {
828: eta_q += Ni[q][k] * eta_field_e[k];
829: rho_q += Ni[q][k] * rho_field_e[k];
830: }
832: q_eta[nqp * e + q] = PetscRealPart(eta_q);
833: q_rhs[nqp * e + q] = PetscRealPart(rho_q);
834: }
835: }
836: PetscCall(DMDARestoreElements(dm, &nel, &npe, &element_list));
838: PetscCall(DMSwarmRestoreField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
839: PetscCall(DMSwarmRestoreField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
841: PetscCall(VecRestoreArray(rho_l, &_rho_l));
842: PetscCall(VecRestoreArray(eta_l, &_eta_l));
843: PetscCall(DMRestoreLocalVector(dm, &rho_l));
844: PetscCall(DMRestoreLocalVector(dm, &eta_l));
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
848: static PetscErrorCode SolveTimeDepStokes(PetscInt mx, PetscInt my)
849: {
850: DM dm_stokes, dm_coeff;
851: PetscInt u_dof, p_dof, dof, stencil_width;
852: Mat A, B;
853: PetscInt nel_local;
854: Vec eta_v, rho_v;
855: Vec f, X;
856: KSP ksp;
857: PC pc;
858: char filename[PETSC_MAX_PATH_LEN];
859: DM dms_quadrature, dms_mpoint;
860: PetscInt nel, npe, npoints;
861: const PetscInt *element_list;
862: PetscInt tk, nt, dump_freq;
863: PetscReal dt, dt_max = 0.0;
864: PetscReal vx[2], vy[2], max_v = 0.0, max_v_step, dh;
865: const char *fieldnames[] = {"eta", "rho"};
866: Vec *pfields;
867: PetscInt ppcell = 1;
868: PetscReal time, delta_eta = 1.0;
869: PetscBool randomize_coords = PETSC_FALSE;
870: PetscReal randomize_fac = 0.25;
871: PetscBool no_view = PETSC_FALSE;
872: PetscBool isbddc;
874: PetscFunctionBeginUser;
875: /*
876: Generate the DMDA for the velocity and pressure spaces.
877: We use Q1 elements for both fields.
878: The Q1 FE basis on a regular mesh has a 9-point stencil (DMDA_STENCIL_BOX)
879: The number of nodes in each direction is mx+1, my+1
880: */
881: u_dof = U_DOFS; /* Vx, Vy - velocities */
882: p_dof = P_DOFS; /* p - pressure */
883: dof = u_dof + p_dof;
884: stencil_width = 1;
885: 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, &dm_stokes));
886: PetscCall(DMDASetElementType(dm_stokes, DMDA_ELEMENT_Q1));
887: PetscCall(DMSetMatType(dm_stokes, MATAIJ));
888: PetscCall(DMSetFromOptions(dm_stokes));
889: PetscCall(DMSetUp(dm_stokes));
890: PetscCall(DMDASetFieldName(dm_stokes, 0, "ux"));
891: PetscCall(DMDASetFieldName(dm_stokes, 1, "uy"));
892: PetscCall(DMDASetFieldName(dm_stokes, 2, "p"));
894: /* unit box [0,0.9142] x [0,1] */
895: PetscCall(DMDASetUniformCoordinates(dm_stokes, 0.0, 0.9142, 0.0, 1.0, 0., 0.));
896: dh = 1.0 / ((PetscReal)(mx));
898: /* Get local number of elements */
899: {
900: PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));
902: nel_local = nel;
904: PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
905: }
907: /* Create DMDA for representing scalar fields */
908: PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 1, &dm_coeff));
910: /* Create the swarm for storing quadrature point values */
911: PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_quadrature));
912: PetscCall(DMSetType(dms_quadrature, DMSWARM));
913: PetscCall(DMSetDimension(dms_quadrature, 2));
914: PetscCall(PetscObjectSetName((PetscObject)dms_quadrature, "Quadrature Swarm"));
916: /* Register fields for viscosity and density on the quadrature points */
917: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "eta_q", 1, PETSC_REAL));
918: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "rho_q", 1, PETSC_REAL));
919: PetscCall(DMSwarmFinalizeFieldRegister(dms_quadrature));
920: PetscCall(DMSwarmSetLocalSizes(dms_quadrature, nel_local * GAUSS_POINTS, 0));
922: /* Create the material point swarm */
923: PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_mpoint));
924: PetscCall(DMSetType(dms_mpoint, DMSWARM));
925: PetscCall(DMSetDimension(dms_mpoint, 2));
926: PetscCall(PetscObjectSetName((PetscObject)dms_mpoint, "Material Point Swarm"));
928: /* Configure the material point swarm to be of type Particle-In-Cell */
929: PetscCall(DMSwarmSetType(dms_mpoint, DMSWARM_PIC));
931: /*
932: Specify the DM to use for point location and projections
933: within the context of a PIC scheme
934: */
935: PetscCall(DMSwarmSetCellDM(dms_mpoint, dm_coeff));
937: /* Register fields for viscosity and density */
938: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "eta", 1, PETSC_REAL));
939: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "rho", 1, PETSC_REAL));
940: PetscCall(DMSwarmFinalizeFieldRegister(dms_mpoint));
942: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL));
943: PetscCall(DMSwarmSetLocalSizes(dms_mpoint, nel_local * ppcell, 100));
945: /*
946: Layout the material points in space using the cell DM.
947: Particle coordinates are defined by cell wise using different methods.
948: - DMSWARMPIC_LAYOUT_GAUSS defines particles coordinates at the positions
949: corresponding to a Gauss quadrature rule with
950: ppcell points in each direction.
951: - DMSWARMPIC_LAYOUT_REGULAR defines particle coordinates at the centoid of
952: ppcell x ppcell quadralaterals defined within the
953: reference element.
954: - DMSWARMPIC_LAYOUT_SUBDIVISION defines particles coordinates at the centroid
955: of each quadralateral obtained by sub-dividing
956: the reference element cell ppcell times.
957: */
958: PetscCall(DMSwarmInsertPointsUsingCellDM(dms_mpoint, DMSWARMPIC_LAYOUT_SUBDIVISION, ppcell));
960: /*
961: Defne a high resolution layer of material points across the material interface
962: */
963: {
964: PetscInt npoints_dir_x[2];
965: PetscReal min[2], max[2];
967: npoints_dir_x[0] = (PetscInt)(0.9142 / (0.05 * dh));
968: npoints_dir_x[1] = (PetscInt)((0.25 - 0.15) / (0.05 * dh));
969: min[0] = 0.0;
970: max[0] = 0.9142;
971: min[1] = 0.05;
972: max[1] = 0.35;
973: PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
974: }
976: /*
977: Define a high resolution layer of material points near the surface of the domain
978: to deal with weakly compressible Q1-Q1 elements. These elements "self compact"
979: when applied to buouyancy driven flow. The error in div(u) is O(h).
980: */
981: {
982: PetscInt npoints_dir_x[2];
983: PetscReal min[2], max[2];
985: npoints_dir_x[0] = (PetscInt)(0.9142 / (0.25 * dh));
986: npoints_dir_x[1] = (PetscInt)(3.0 * dh / (0.25 * dh));
987: min[0] = 0.0;
988: max[0] = 0.9142;
989: min[1] = 1.0 - 3.0 * dh;
990: max[1] = 1.0 - 0.0001;
991: PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
992: }
994: PetscCall(DMView(dms_mpoint, PETSC_VIEWER_STDOUT_WORLD));
996: /* Define initial material properties on each particle in the material point swarm */
997: PetscCall(PetscOptionsGetReal(NULL, NULL, "-delta_eta", &delta_eta, NULL));
998: PetscCall(PetscOptionsGetBool(NULL, NULL, "-randomize_coords", &randomize_coords, NULL));
999: PetscCall(PetscOptionsGetReal(NULL, NULL, "-randomize_fac", &randomize_fac, NULL));
1000: PetscCheck(randomize_fac <= 1.0, PETSC_COMM_WORLD, PETSC_ERR_USER, "The value of -randomize_fac should be <= 1.0");
1001: {
1002: PetscReal *array_x, *array_e, *array_r;
1003: PetscInt p;
1004: PetscRandom r;
1005: PetscMPIInt rank;
1007: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1009: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
1010: PetscCall(PetscRandomSetInterval(r, -randomize_fac * dh, randomize_fac * dh));
1011: PetscCall(PetscRandomSetSeed(r, (unsigned long)rank));
1012: PetscCall(PetscRandomSeed(r));
1014: PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));
1016: /*
1017: Fetch the registered data from the material point DMSwarm.
1018: The fields "eta" and "rho" were registered by this example.
1019: The field identified by the the variable DMSwarmPICField_coor
1020: was registered by the DMSwarm implementation when the function
1021: DMSwarmSetType(dms_mpoint,DMSWARM_PIC)
1022: was called. The returned array defines the coordinates of each
1023: material point in the point swarm.
1024: */
1025: PetscCall(DMSwarmGetField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));
1026: PetscCall(DMSwarmGetField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
1027: PetscCall(DMSwarmGetField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));
1029: PetscCall(DMSwarmGetLocalSize(dms_mpoint, &npoints));
1030: for (p = 0; p < npoints; p++) {
1031: PetscReal x_p[2], rr[2];
1033: if (randomize_coords) {
1034: PetscCall(PetscRandomGetValueReal(r, &rr[0]));
1035: PetscCall(PetscRandomGetValueReal(r, &rr[1]));
1036: array_x[2 * p + 0] += rr[0];
1037: array_x[2 * p + 1] += rr[1];
1038: }
1040: /* Get the coordinates of point, p */
1041: x_p[0] = array_x[2 * p + 0];
1042: x_p[1] = array_x[2 * p + 1];
1044: if (x_p[1] < (0.2 + 0.02 * PetscCosReal(PETSC_PI * x_p[0] / 0.9142))) {
1045: /* Material properties below the interface */
1046: array_e[p] = 1.0 * (1.0 / delta_eta);
1047: array_r[p] = 0.0;
1048: } else {
1049: /* Material properties above the interface */
1050: array_e[p] = 1.0;
1051: array_r[p] = 1.0;
1052: }
1053: }
1055: /*
1056: Restore the fetched data fields from the material point DMSwarm.
1057: Calling the Restore function invalidates the points array_r, array_e, array_x
1058: by setting them to NULL.
1059: */
1060: PetscCall(DMSwarmRestoreField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));
1061: PetscCall(DMSwarmRestoreField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
1062: PetscCall(DMSwarmRestoreField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));
1064: PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
1065: PetscCall(PetscRandomDestroy(&r));
1066: }
1068: /*
1069: If the particle coordinates where randomly shifted, they may have crossed into another
1070: element, or into another sub-domain. To account for this we call the Migrate function.
1071: */
1072: if (randomize_coords) PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));
1074: PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL));
1075: if (!no_view) PetscCall(DMSwarmViewXDMF(dms_mpoint, "ic_coeff_dms.xmf"));
1077: /* project the swarm properties */
1078: PetscCall(DMSwarmProjectFields(dms_mpoint, 2, fieldnames, &pfields, PETSC_FALSE));
1079: eta_v = pfields[0];
1080: rho_v = pfields[1];
1081: PetscCall(PetscObjectSetName((PetscObject)eta_v, "eta"));
1082: PetscCall(PetscObjectSetName((PetscObject)rho_v, "rho"));
1083: PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));
1085: /* view projected coefficients eta and rho */
1086: if (!no_view) {
1087: PetscViewer viewer;
1089: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1090: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1091: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1092: PetscCall(PetscViewerFileSetName(viewer, "ic_coeff_dmda.vts"));
1093: PetscCall(VecView(eta_v, viewer));
1094: PetscCall(VecView(rho_v, viewer));
1095: PetscCall(PetscViewerDestroy(&viewer));
1096: }
1098: PetscCall(DMCreateMatrix(dm_stokes, &A));
1099: PetscCall(DMCreateMatrix(dm_stokes, &B));
1100: PetscCall(DMCreateGlobalVector(dm_stokes, &f));
1101: PetscCall(DMCreateGlobalVector(dm_stokes, &X));
1103: PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
1104: PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
1105: PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));
1107: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
1108: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));
1110: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
1111: PetscCall(KSPSetOptionsPrefix(ksp, "stokes_"));
1112: PetscCall(KSPSetDM(ksp, dm_stokes));
1113: PetscCall(KSPSetDMActive(ksp, PETSC_FALSE));
1114: PetscCall(KSPSetOperators(ksp, A, B));
1115: PetscCall(KSPSetFromOptions(ksp));
1116: PetscCall(KSPGetPC(ksp, &pc));
1117: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &isbddc));
1118: if (isbddc) PetscCall(KSPSetOperators(ksp, A, A));
1120: /* Define u-v-p indices for fieldsplit */
1121: {
1122: PC pc;
1123: const PetscInt ufields[] = {0, 1}, pfields[1] = {2};
1125: PetscCall(KSPGetPC(ksp, &pc));
1126: PetscCall(PCFieldSplitSetBlockSize(pc, 3));
1127: PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
1128: PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
1129: }
1131: /* If using a fieldsplit preconditioner, attach a DMDA to the velocity split so that geometric multigrid can be used */
1132: {
1133: PC pc, pc_u;
1134: KSP *sub_ksp, ksp_u;
1135: PetscInt nsplits;
1136: DM dm_u;
1137: PetscBool is_pcfs;
1139: PetscCall(KSPGetPC(ksp, &pc));
1141: is_pcfs = PETSC_FALSE;
1142: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs));
1144: if (is_pcfs) {
1145: PetscCall(KSPSetUp(ksp));
1146: PetscCall(KSPGetPC(ksp, &pc));
1147: PetscCall(PCFieldSplitGetSubKSP(pc, &nsplits, &sub_ksp));
1148: ksp_u = sub_ksp[0];
1149: PetscCall(PetscFree(sub_ksp));
1151: if (nsplits == 2) {
1152: PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 2, &dm_u));
1154: PetscCall(KSPSetDM(ksp_u, dm_u));
1155: PetscCall(KSPSetDMActive(ksp_u, PETSC_FALSE));
1156: PetscCall(DMDestroy(&dm_u));
1158: /* enforce galerkin coarse grids be used */
1159: PetscCall(KSPGetPC(ksp_u, &pc_u));
1160: PetscCall(PCMGSetGalerkin(pc_u, PC_MG_GALERKIN_PMAT));
1161: }
1162: }
1163: }
1165: dump_freq = 10;
1166: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dump_freq", &dump_freq, NULL));
1167: nt = 10;
1168: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nt", &nt, NULL));
1169: time = 0.0;
1170: for (tk = 1; tk <= nt; tk++) {
1171: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... assemble\n"));
1172: PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
1173: PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
1174: PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));
1176: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... bc imposition\n"));
1177: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
1178: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));
1180: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... solve\n"));
1181: PetscCall(KSPSetOperators(ksp, A, isbddc ? A : B));
1182: PetscCall(KSPSolve(ksp, f, X));
1184: PetscCall(VecStrideMax(X, 0, NULL, &vx[1]));
1185: PetscCall(VecStrideMax(X, 1, NULL, &vy[1]));
1186: PetscCall(VecStrideMin(X, 0, NULL, &vx[0]));
1187: PetscCall(VecStrideMin(X, 1, NULL, &vy[0]));
1189: max_v_step = PetscMax(vx[0], vx[1]);
1190: max_v_step = PetscMax(max_v_step, vy[0]);
1191: max_v_step = PetscMax(max_v_step, vy[1]);
1192: max_v = PetscMax(max_v, max_v_step);
1194: dt_max = 2.0;
1195: dt = 0.5 * (dh / max_v_step);
1196: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... max v %1.4e , dt %1.4e : [total] max v %1.4e , dt_max %1.4e\n", (double)max_v_step, (double)dt, (double)max_v, (double)dt_max));
1197: dt = PetscMin(dt_max, dt);
1199: /* advect */
1200: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... advect\n"));
1201: PetscCall(MaterialPoint_AdvectRK1(dm_stokes, X, dt, dms_mpoint));
1203: /* migrate */
1204: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... migrate\n"));
1205: PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));
1207: /* update cell population */
1208: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... populate cells\n"));
1209: PetscCall(MaterialPoint_PopulateCell(dm_stokes, dms_mpoint));
1211: /* update coefficients on quadrature points */
1212: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... project\n"));
1213: PetscCall(DMSwarmProjectFields(dms_mpoint, 2, fieldnames, &pfields, PETSC_TRUE));
1214: eta_v = pfields[0];
1215: rho_v = pfields[1];
1216: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... interp\n"));
1217: PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));
1219: if (tk % dump_freq == 0) {
1220: PetscViewer viewer;
1222: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... write XDMF, VTS\n"));
1223: PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_coeff_dms.xmf", tk));
1224: PetscCall(DMSwarmViewXDMF(dms_mpoint, filename));
1226: PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_vp_dm.vts", tk));
1227: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1228: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1229: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1230: PetscCall(PetscViewerFileSetName(viewer, filename));
1231: PetscCall(VecView(X, viewer));
1232: PetscCall(PetscViewerDestroy(&viewer));
1233: }
1234: time += dt;
1235: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step %" PetscInt_FMT " : time %1.2e \n", tk, (double)time));
1236: }
1238: PetscCall(KSPDestroy(&ksp));
1239: PetscCall(VecDestroy(&X));
1240: PetscCall(VecDestroy(&f));
1241: PetscCall(MatDestroy(&A));
1242: PetscCall(MatDestroy(&B));
1243: PetscCall(VecDestroy(&eta_v));
1244: PetscCall(VecDestroy(&rho_v));
1245: PetscCall(PetscFree(pfields));
1247: PetscCall(DMDestroy(&dms_mpoint));
1248: PetscCall(DMDestroy(&dms_quadrature));
1249: PetscCall(DMDestroy(&dm_coeff));
1250: PetscCall(DMDestroy(&dm_stokes));
1251: PetscFunctionReturn(PETSC_SUCCESS);
1252: }
1254: /*
1255: <sequential run>
1256: ./ex70 -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 lu -mx 80 -my 80 -stokes_ksp_converged_reason -dump_freq 25 -stokes_ksp_rtol 1.0e-8 -build_twosided allreduce -ppcell 2 -nt 4000 -delta_eta 1.0 -randomize_coords
1257: */
1258: int main(int argc, char **args)
1259: {
1260: PetscInt mx, my;
1261: PetscBool set = PETSC_FALSE;
1263: PetscFunctionBeginUser;
1264: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
1265: mx = my = 10;
1266: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1267: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1268: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mxy", &mx, &set));
1269: if (set) my = mx;
1270: PetscCall(SolveTimeDepStokes(mx, my));
1271: PetscCall(PetscFinalize());
1272: return 0;
1273: }
1275: /* -------------------------- helpers for boundary conditions -------------------------------- */
1276: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
1277: {
1278: DM cda;
1279: Vec coords;
1280: PetscInt si, sj, nx, ny, i, j;
1281: PetscInt M, N;
1282: DMDACoor2d **_coords;
1283: const PetscInt *g_idx;
1284: PetscInt *bc_global_ids;
1285: PetscScalar *bc_vals;
1286: PetscInt nbcs;
1287: PetscInt n_dofs;
1288: ISLocalToGlobalMapping ltogm;
1290: PetscFunctionBeginUser;
1291: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1292: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1294: PetscCall(DMGetCoordinateDM(da, &cda));
1295: PetscCall(DMGetCoordinatesLocal(da, &coords));
1296: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1297: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1298: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1300: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1301: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1303: /* init the entries to -1 so VecSetValues will ignore them */
1304: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1306: i = nx - 1;
1307: for (j = 0; j < ny; j++) {
1308: PetscInt local_id;
1310: local_id = i + j * nx;
1312: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1314: bc_vals[j] = 0.0;
1315: }
1316: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1317: nbcs = 0;
1318: if ((si + nx) == (M)) nbcs = ny;
1320: if (b) {
1321: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1322: PetscCall(VecAssemblyBegin(b));
1323: PetscCall(VecAssemblyEnd(b));
1324: }
1325: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1327: PetscCall(PetscFree(bc_vals));
1328: PetscCall(PetscFree(bc_global_ids));
1330: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1331: PetscFunctionReturn(PETSC_SUCCESS);
1332: }
1334: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
1335: {
1336: DM cda;
1337: Vec coords;
1338: PetscInt si, sj, nx, ny, i, j;
1339: PetscInt M, N;
1340: DMDACoor2d **_coords;
1341: const PetscInt *g_idx;
1342: PetscInt *bc_global_ids;
1343: PetscScalar *bc_vals;
1344: PetscInt nbcs;
1345: PetscInt n_dofs;
1346: ISLocalToGlobalMapping ltogm;
1348: PetscFunctionBeginUser;
1349: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1350: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1352: PetscCall(DMGetCoordinateDM(da, &cda));
1353: PetscCall(DMGetCoordinatesLocal(da, &coords));
1354: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1355: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1356: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1358: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1359: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1361: /* init the entries to -1 so VecSetValues will ignore them */
1362: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1364: i = 0;
1365: for (j = 0; j < ny; j++) {
1366: PetscInt local_id;
1368: local_id = i + j * nx;
1370: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1372: bc_vals[j] = 0.0;
1373: }
1374: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1375: nbcs = 0;
1376: if (si == 0) nbcs = ny;
1378: if (b) {
1379: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1380: PetscCall(VecAssemblyBegin(b));
1381: PetscCall(VecAssemblyEnd(b));
1382: }
1384: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1386: PetscCall(PetscFree(bc_vals));
1387: PetscCall(PetscFree(bc_global_ids));
1389: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1390: PetscFunctionReturn(PETSC_SUCCESS);
1391: }
1393: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
1394: {
1395: DM cda;
1396: Vec coords;
1397: PetscInt si, sj, nx, ny, i, j;
1398: PetscInt M, N;
1399: DMDACoor2d **_coords;
1400: const PetscInt *g_idx;
1401: PetscInt *bc_global_ids;
1402: PetscScalar *bc_vals;
1403: PetscInt nbcs;
1404: PetscInt n_dofs;
1405: ISLocalToGlobalMapping ltogm;
1407: PetscFunctionBeginUser;
1408: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1409: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1411: PetscCall(DMGetCoordinateDM(da, &cda));
1412: PetscCall(DMGetCoordinatesLocal(da, &coords));
1413: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1414: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1415: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1417: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1418: PetscCall(PetscMalloc1(nx, &bc_vals));
1420: /* init the entries to -1 so VecSetValues will ignore them */
1421: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1423: j = ny - 1;
1424: for (i = 0; i < nx; i++) {
1425: PetscInt local_id;
1427: local_id = i + j * nx;
1429: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1431: bc_vals[i] = 0.0;
1432: }
1433: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1434: nbcs = 0;
1435: if ((sj + ny) == (N)) nbcs = nx;
1437: if (b) {
1438: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1439: PetscCall(VecAssemblyBegin(b));
1440: PetscCall(VecAssemblyEnd(b));
1441: }
1442: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, NULL, NULL));
1444: PetscCall(PetscFree(bc_vals));
1445: PetscCall(PetscFree(bc_global_ids));
1447: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1448: PetscFunctionReturn(PETSC_SUCCESS);
1449: }
1451: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
1452: {
1453: DM cda;
1454: Vec coords;
1455: PetscInt si, sj, nx, ny, i, j;
1456: PetscInt M, N;
1457: DMDACoor2d **_coords;
1458: const PetscInt *g_idx;
1459: PetscInt *bc_global_ids;
1460: PetscScalar *bc_vals;
1461: PetscInt nbcs;
1462: PetscInt n_dofs;
1463: ISLocalToGlobalMapping ltogm;
1465: PetscFunctionBeginUser;
1466: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1467: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1469: PetscCall(DMGetCoordinateDM(da, &cda));
1470: PetscCall(DMGetCoordinatesLocal(da, &coords));
1471: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1472: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1473: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1475: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1476: PetscCall(PetscMalloc1(nx, &bc_vals));
1478: /* init the entries to -1 so VecSetValues will ignore them */
1479: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1481: j = 0;
1482: for (i = 0; i < nx; i++) {
1483: PetscInt local_id;
1485: local_id = i + j * nx;
1487: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1489: bc_vals[i] = 0.0;
1490: }
1491: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1492: nbcs = 0;
1493: if (sj == 0) nbcs = nx;
1495: if (b) {
1496: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1497: PetscCall(VecAssemblyBegin(b));
1498: PetscCall(VecAssemblyEnd(b));
1499: }
1500: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1502: PetscCall(PetscFree(bc_vals));
1503: PetscCall(PetscFree(bc_global_ids));
1505: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1506: PetscFunctionReturn(PETSC_SUCCESS);
1507: }
1509: /*
1510: Impose free slip boundary conditions on the left/right faces: u_i n_i = 0, tau_{ij} t_j = 0
1511: Impose no slip boundray conditions on the top/bottom faces: u_i n_i = 0, u_i t_i = 0
1512: */
1513: static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes, Mat A, Vec f)
1514: {
1515: PetscFunctionBeginUser;
1516: PetscCall(BCApplyZero_NORTH(dm_stokes, 0, A, f));
1517: PetscCall(BCApplyZero_NORTH(dm_stokes, 1, A, f));
1518: PetscCall(BCApplyZero_EAST(dm_stokes, 0, A, f));
1519: PetscCall(BCApplyZero_SOUTH(dm_stokes, 0, A, f));
1520: PetscCall(BCApplyZero_SOUTH(dm_stokes, 1, A, f));
1521: PetscCall(BCApplyZero_WEST(dm_stokes, 0, A, f));
1522: PetscFunctionReturn(PETSC_SUCCESS);
1523: }
1525: /*TEST
1527: test:
1528: suffix: 1
1529: args: -no_view
1530: requires: !complex double
1531: filter: grep -v atomic
1532: filter_output: grep -v atomic
1533: test:
1534: suffix: 1_matis
1535: requires: !complex double
1536: args: -no_view -dm_mat_type is
1537: filter: grep -v atomic
1538: filter_output: grep -v atomic
1539: testset:
1540: nsize: 4
1541: requires: !complex double
1542: args: -no_view -dm_mat_type is -stokes_ksp_type fetidp -mx 80 -my 80 -stokes_ksp_converged_reason -stokes_ksp_rtol 1.0e-8 -ppcell 2 -nt 4 -randomize_coords -stokes_ksp_error_if_not_converged
1543: filter: grep -v atomic
1544: filter_output: grep -v atomic
1545: test:
1546: suffix: fetidp
1547: args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
1548: test:
1549: suffix: fetidp_lumped
1550: args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd -stokes_fetidp_pc_lumped -stokes_fetidp_bddc_pc_bddc_dirichlet_pc_type none -stokes_fetidp_bddc_pc_bddc_switch_static
1551: test:
1552: suffix: fetidp_saddlepoint
1553: args: -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_norm_type natural -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
1554: test:
1555: suffix: fetidp_saddlepoint_lumped
1556: args: -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_norm_type natural -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 -stokes_fetidp_bddc_pc_bddc_dirichlet_pc_type none -stokes_fetidp_bddc_pc_bddc_switch_static -stokes_fetidp_pc_lumped
1557: TEST*/