Actual source code: ex70.c
1: static char help[] = "------------------------------------------------------------------------------------------------------------------------------ \n\
2: Solves the time-dependent incompressible, variable viscosity Stokes equation in 2D driven by buoyancy 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 viscosity 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>
35: static PetscErrorCode DMDAApplyBoundaryConditions(DM, Mat, Vec);
37: #define NSD 2 /* number of spatial dimensions */
38: #define NODES_PER_EL 4 /* nodes per element */
39: #define U_DOFS 2 /* degrees of freedom per velocity node */
40: #define P_DOFS 1 /* degrees of freedom per pressure node */
41: #define GAUSS_POINTS 4
43: static void EvaluateBasis_Q1(PetscScalar _xi[], PetscScalar N[])
44: {
45: PetscScalar xi = _xi[0];
46: PetscScalar eta = _xi[1];
48: N[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
49: N[1] = 0.25 * (1.0 + xi) * (1.0 - eta);
50: N[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
51: N[3] = 0.25 * (1.0 - xi) * (1.0 + eta);
52: }
54: static void EvaluateBasisDerivatives_Q1(PetscScalar _xi[], PetscScalar dN[][NODES_PER_EL])
55: {
56: PetscScalar xi = _xi[0];
57: PetscScalar eta = _xi[1];
59: dN[0][0] = -0.25 * (1.0 - eta);
60: dN[0][1] = 0.25 * (1.0 - eta);
61: dN[0][2] = 0.25 * (1.0 + eta);
62: dN[0][3] = -0.25 * (1.0 + eta);
64: dN[1][0] = -0.25 * (1.0 - xi);
65: dN[1][1] = -0.25 * (1.0 + xi);
66: dN[1][2] = 0.25 * (1.0 + xi);
67: dN[1][3] = 0.25 * (1.0 - xi);
68: }
70: static void EvaluateDerivatives(PetscScalar dN[][NODES_PER_EL], PetscScalar dNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
71: {
72: PetscScalar J00, J01, J10, J11, J;
73: PetscScalar iJ00, iJ01, iJ10, iJ11;
74: PetscInt i;
76: J00 = J01 = J10 = J11 = 0.0;
77: for (i = 0; i < NODES_PER_EL; i++) {
78: PetscScalar cx = coords[2 * i];
79: PetscScalar cy = coords[2 * i + 1];
81: J00 += dN[0][i] * cx; /* J_xx = dx/dxi */
82: J01 += dN[0][i] * cy; /* J_xy = dy/dxi */
83: J10 += dN[1][i] * cx; /* J_yx = dx/deta */
84: J11 += dN[1][i] * cy; /* J_yy = dy/deta */
85: }
86: J = (J00 * J11) - (J01 * J10);
88: iJ00 = J11 / J;
89: iJ01 = -J01 / J;
90: iJ10 = -J10 / J;
91: iJ11 = J00 / J;
93: for (i = 0; i < NODES_PER_EL; i++) {
94: dNx[0][i] = dN[0][i] * iJ00 + dN[1][i] * iJ01;
95: dNx[1][i] = dN[0][i] * iJ10 + dN[1][i] * iJ11;
96: }
98: if (det_J) *det_J = J;
99: }
101: static void CreateGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
102: {
103: *ngp = 4;
104: gp_xi[0][0] = -0.57735026919;
105: gp_xi[0][1] = -0.57735026919;
106: gp_xi[1][0] = -0.57735026919;
107: gp_xi[1][1] = 0.57735026919;
108: gp_xi[2][0] = 0.57735026919;
109: gp_xi[2][1] = 0.57735026919;
110: gp_xi[3][0] = 0.57735026919;
111: gp_xi[3][1] = -0.57735026919;
112: gp_weight[0] = 1.0;
113: gp_weight[1] = 1.0;
114: gp_weight[2] = 1.0;
115: gp_weight[3] = 1.0;
116: }
118: static PetscErrorCode DMDAGetElementEqnums_up(const PetscInt element[], PetscInt s_u[], PetscInt s_p[])
119: {
120: PetscFunctionBeginUser;
121: for (PetscInt i = 0; i < NODES_PER_EL; i++) {
122: /* velocity */
123: s_u[NSD * i + 0] = 3 * element[i];
124: s_u[NSD * i + 1] = 3 * element[i] + 1;
125: /* pressure */
126: s_p[i] = 3 * element[i] + 2;
127: }
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: 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)
132: {
133: PetscInt ij, r, c, nc;
135: nc = u_NPE * u_dof;
136: r = w_dof * wi + wd;
137: c = u_dof * ui + ud;
138: ij = r * nc + c;
139: return ij;
140: }
142: static void BForm_DivT(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
143: {
144: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
145: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
146: PetscScalar J_p, tildeD[3];
147: PetscScalar B[3][U_DOFS * NODES_PER_EL];
148: PetscInt p, i, j, k, ngp;
150: /* define quadrature rule */
151: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
153: /* evaluate bilinear form */
154: for (p = 0; p < ngp; p++) {
155: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
156: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
158: for (i = 0; i < NODES_PER_EL; i++) {
159: PetscScalar d_dx_i = GNx_p[0][i];
160: PetscScalar d_dy_i = GNx_p[1][i];
162: B[0][2 * i] = d_dx_i;
163: B[0][2 * i + 1] = 0.0;
164: B[1][2 * i] = 0.0;
165: B[1][2 * i + 1] = d_dy_i;
166: B[2][2 * i] = d_dy_i;
167: B[2][2 * i + 1] = d_dx_i;
168: }
170: tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
171: tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
172: tildeD[2] = gp_weight[p] * J_p * eta[p];
174: /* form Bt tildeD B */
175: /*
176: Ke_ij = Bt_ik . D_kl . B_lj
177: = B_ki . D_kl . B_lj
178: = B_ki . D_kk . B_kj
179: */
180: for (i = 0; i < 8; i++) {
181: for (j = 0; j < 8; j++) {
182: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
183: Ke[i + 8 * j] += B[k][i] * tildeD[k] * B[k][j];
184: }
185: }
186: }
187: }
188: }
190: static void BForm_Grad(PetscScalar Ke[], PetscScalar coords[])
191: {
192: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
193: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
194: PetscScalar J_p, fac;
195: PetscInt p, i, j, di, ngp;
197: /* define quadrature rule */
198: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
200: /* evaluate bilinear form */
201: for (p = 0; p < ngp; p++) {
202: EvaluateBasis_Q1(gp_xi[p], Ni_p);
203: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
204: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
205: fac = gp_weight[p] * J_p;
207: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
208: for (di = 0; di < NSD; di++) { /* u dofs */
209: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
210: PetscInt IJ;
211: IJ = map_wIwDI_uJuDJ(i, di, NODES_PER_EL, 2, j, 0, NODES_PER_EL, 1);
213: Ke[IJ] -= GNx_p[di][i] * Ni_p[j] * fac;
214: }
215: }
216: }
217: }
218: }
220: static void BForm_Div(PetscScalar De[], PetscScalar coords[])
221: {
222: PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
223: PetscInt i, j, nr_g, nc_g;
225: PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));
226: BForm_Grad(Ge, coords);
228: nr_g = U_DOFS * NODES_PER_EL;
229: nc_g = P_DOFS * NODES_PER_EL;
231: for (i = 0; i < nr_g; i++) {
232: for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
233: }
234: }
236: static void BForm_Stabilisation(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
237: {
238: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
239: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
240: PetscScalar J_p, fac, eta_avg;
241: PetscInt p, i, j, ngp;
243: /* define quadrature rule */
244: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
246: /* evaluate bilinear form */
247: for (p = 0; p < ngp; p++) {
248: EvaluateBasis_Q1(gp_xi[p], Ni_p);
249: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
250: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
251: fac = gp_weight[p] * J_p;
253: for (i = 0; i < NODES_PER_EL; i++) {
254: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * (Ni_p[i] * Ni_p[j] - 0.0625);
255: }
256: }
258: /* scale */
259: eta_avg = 0.0;
260: for (p = 0; p < ngp; p++) eta_avg += eta[p];
261: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
262: fac = 1.0 / eta_avg;
263: for (i = 0; i < NODES_PER_EL; i++) {
264: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
265: }
266: }
268: static void BForm_ScaledMassMatrix(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
269: {
270: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
271: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
272: PetscScalar J_p, fac, eta_avg;
273: PetscInt p, i, j, ngp;
275: /* define quadrature rule */
276: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
278: /* evaluate bilinear form */
279: for (p = 0; p < ngp; p++) {
280: EvaluateBasis_Q1(gp_xi[p], Ni_p);
281: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
282: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
283: fac = gp_weight[p] * J_p;
285: for (i = 0; i < NODES_PER_EL; i++) {
286: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * Ni_p[i] * Ni_p[j];
287: }
288: }
290: /* scale */
291: eta_avg = 0.0;
292: for (p = 0; p < ngp; p++) eta_avg += eta[p];
293: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
294: fac = 1.0 / eta_avg;
295: for (i = 0; i < NODES_PER_EL; i++) {
296: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] *= fac;
297: }
298: }
300: static void LForm_MomentumRHS(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
301: {
302: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
303: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
304: PetscScalar J_p, fac;
305: PetscInt p, i, ngp;
307: /* define quadrature rule */
308: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
310: /* evaluate linear form */
311: for (p = 0; p < ngp; p++) {
312: EvaluateBasis_Q1(gp_xi[p], Ni_p);
313: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
314: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
315: fac = gp_weight[p] * J_p;
317: for (i = 0; i < NODES_PER_EL; i++) {
318: Fe[NSD * i] = 0.0;
319: Fe[NSD * i + 1] -= fac * Ni_p[i] * fy[p];
320: }
321: }
322: }
324: static PetscErrorCode GetElementCoords(const PetscScalar _coords[], const PetscInt e2n[], PetscScalar el_coords[])
325: {
326: PetscFunctionBeginUser;
327: /* get coords for the element */
328: for (PetscInt i = 0; i < 4; i++) {
329: for (PetscInt d = 0; d < NSD; d++) el_coords[NSD * i + d] = _coords[NSD * e2n[i] + d];
330: }
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: static PetscErrorCode AssembleStokes_A(Mat A, DM stokes_da, DM quadrature)
335: {
336: DM cda;
337: Vec coords;
338: const PetscScalar *_coords;
339: PetscInt u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
340: PetscInt p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
341: PetscInt nel, npe, eidx;
342: const PetscInt *element_list;
343: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
344: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
345: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
346: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
347: PetscScalar el_coords[NODES_PER_EL * NSD];
348: PetscScalar *q_eta, *prop_eta;
350: PetscFunctionBeginUser;
351: PetscCall(MatZeroEntries(A));
352: /* setup for coords */
353: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
354: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
355: PetscCall(VecGetArrayRead(coords, &_coords));
357: /* setup for coefficients */
358: PetscCall(DMSwarmGetField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
360: PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
361: for (eidx = 0; eidx < nel; eidx++) {
362: const PetscInt *element = &element_list[npe * eidx];
364: /* get coords for the element */
365: PetscCall(GetElementCoords(_coords, element, el_coords));
367: /* get coefficients for the element */
368: prop_eta = &q_eta[GAUSS_POINTS * eidx];
370: /* initialise element stiffness matrix */
371: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
372: PetscCall(PetscMemzero(Ge, sizeof(Ge)));
373: PetscCall(PetscMemzero(De, sizeof(De)));
374: PetscCall(PetscMemzero(Ce, sizeof(Ce)));
376: /* form element stiffness matrix */
377: BForm_DivT(Ae, el_coords, prop_eta);
378: BForm_Grad(Ge, el_coords);
379: BForm_Div(De, el_coords);
380: BForm_Stabilisation(Ce, el_coords, prop_eta);
382: /* insert element matrix into global matrix */
383: PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
384: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
385: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
386: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
387: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
388: }
389: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
390: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
392: PetscCall(DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
393: PetscCall(VecRestoreArrayRead(coords, &_coords));
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: static PetscErrorCode AssembleStokes_PC(Mat A, DM stokes_da, DM quadrature)
398: {
399: DM cda;
400: Vec coords;
401: const PetscScalar *_coords;
402: PetscInt u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
403: PetscInt p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
404: PetscInt nel, npe, eidx;
405: const PetscInt *element_list;
406: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
407: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
408: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
409: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
410: PetscScalar el_coords[NODES_PER_EL * NSD];
411: PetscScalar *q_eta, *prop_eta;
413: PetscFunctionBeginUser;
414: PetscCall(MatZeroEntries(A));
415: /* setup for coords */
416: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
417: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
418: PetscCall(VecGetArrayRead(coords, &_coords));
420: /* setup for coefficients */
421: PetscCall(DMSwarmGetField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
423: PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
424: for (eidx = 0; eidx < nel; eidx++) {
425: const PetscInt *element = &element_list[npe * eidx];
427: /* get coords for the element */
428: PetscCall(GetElementCoords(_coords, element, el_coords));
430: /* get coefficients for the element */
431: prop_eta = &q_eta[GAUSS_POINTS * eidx];
433: /* initialise element stiffness matrix */
434: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
435: PetscCall(PetscMemzero(Ge, sizeof(Ge)));
436: PetscCall(PetscMemzero(De, sizeof(De)));
437: PetscCall(PetscMemzero(Ce, sizeof(Ce)));
439: /* form element stiffness matrix */
440: BForm_DivT(Ae, el_coords, prop_eta);
441: BForm_Grad(Ge, el_coords);
442: BForm_ScaledMassMatrix(Ce, el_coords, prop_eta);
444: /* insert element matrix into global matrix */
445: PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
446: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
447: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
448: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
449: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
450: }
451: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
452: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
454: PetscCall(DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
455: PetscCall(VecRestoreArrayRead(coords, &_coords));
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: static PetscErrorCode AssembleStokes_RHS(Vec F, DM stokes_da, DM quadrature)
460: {
461: DM cda;
462: Vec coords;
463: const PetscScalar *_coords;
464: PetscInt u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
465: PetscInt p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
466: PetscInt nel, npe, eidx, i;
467: const PetscInt *element_list;
468: PetscScalar Fe[NODES_PER_EL * U_DOFS];
469: PetscScalar He[NODES_PER_EL * P_DOFS];
470: PetscScalar el_coords[NODES_PER_EL * NSD];
471: PetscScalar *q_rhs, *prop_fy;
472: Vec local_F;
473: PetscScalar *LA_F;
475: PetscFunctionBeginUser;
476: PetscCall(VecZeroEntries(F));
477: /* setup for coords */
478: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
479: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
480: PetscCall(VecGetArrayRead(coords, &_coords));
482: /* setup for coefficients */
483: PetscCall(DMSwarmGetField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
485: /* get access to the vector */
486: PetscCall(DMGetLocalVector(stokes_da, &local_F));
487: PetscCall(VecZeroEntries(local_F));
488: PetscCall(VecGetArray(local_F, &LA_F));
490: PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
491: for (eidx = 0; eidx < nel; eidx++) {
492: const PetscInt *element = &element_list[npe * eidx];
494: /* get coords for the element */
495: PetscCall(GetElementCoords(_coords, element, el_coords));
497: /* get coefficients for the element */
498: prop_fy = &q_rhs[GAUSS_POINTS * eidx];
500: /* initialise element stiffness matrix */
501: PetscCall(PetscMemzero(Fe, sizeof(Fe)));
502: PetscCall(PetscMemzero(He, sizeof(He)));
504: /* form element stiffness matrix */
505: LForm_MomentumRHS(Fe, el_coords, NULL, prop_fy);
507: /* insert element matrix into global matrix */
508: PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
510: for (i = 0; i < NODES_PER_EL * U_DOFS; i++) LA_F[u_eqn[i]] += Fe[i];
511: }
512: PetscCall(DMSwarmRestoreField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
513: PetscCall(VecRestoreArrayRead(coords, &_coords));
515: PetscCall(VecRestoreArray(local_F, &LA_F));
516: PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
517: PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
518: PetscCall(DMRestoreLocalVector(stokes_da, &local_F));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: PetscErrorCode DMSwarmPICInsertPointsCellwise(DM dm, DM dmc, PetscInt e, PetscInt npoints, PetscReal xi[], PetscBool proximity_initialization)
523: {
524: DMSwarmCellDM celldm;
525: PetscInt dim, nel, npe, q, k, d, ncurr, Nfc;
526: const PetscInt *element_list;
527: Vec coor;
528: const PetscScalar *_coor;
529: PetscReal **basis, *elcoor, *xp;
530: PetscReal *swarm_coor;
531: PetscInt *swarm_cellid;
532: const char **coordFields, *cellid;
534: PetscFunctionBeginUser;
535: PetscCall(DMGetDimension(dm, &dim));
536: PetscCall(DMDAGetElements(dmc, &nel, &npe, &element_list));
538: PetscCall(PetscMalloc1(dim * npoints, &xp));
539: PetscCall(PetscMalloc1(dim * npe, &elcoor));
540: PetscCall(PetscMalloc1(npoints, &basis));
541: for (q = 0; q < npoints; q++) {
542: PetscCall(PetscMalloc1(npe, &basis[q]));
544: switch (dim) {
545: case 1:
546: basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]);
547: basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]);
548: break;
549: case 2:
550: basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
551: basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
552: basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
553: basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
554: break;
556: case 3:
557: basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
558: basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
559: basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
560: basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
561: basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
562: basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
563: basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
564: basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
565: break;
566: }
567: }
569: PetscCall(DMGetCoordinatesLocal(dmc, &coor));
570: PetscCall(VecGetArrayRead(coor, &_coor));
571: /* compute and store the coordinates for the new points */
572: {
573: const PetscInt *element = &element_list[npe * e];
575: for (k = 0; k < npe; k++) {
576: for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]);
577: }
578: for (q = 0; q < npoints; q++) {
579: for (d = 0; d < dim; d++) xp[dim * q + d] = 0.0;
580: for (k = 0; k < npe; k++) {
581: for (d = 0; d < dim; d++) xp[dim * q + d] += basis[q][k] * elcoor[dim * k + d];
582: }
583: }
584: }
585: PetscCall(VecRestoreArrayRead(coor, &_coor));
586: PetscCall(DMDARestoreElements(dmc, &nel, &npe, &element_list));
588: PetscCall(DMSwarmGetLocalSize(dm, &ncurr));
589: PetscCall(DMSwarmAddNPoints(dm, npoints));
590: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
591: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
592: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
593: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
595: if (proximity_initialization) {
596: PetscInt *nnlist;
597: PetscReal *coor_q, *coor_qn;
598: PetscInt npoints_e, *plist_e;
600: PetscCall(DMSwarmSortGetPointsPerCell(dm, e, &npoints_e, &plist_e));
602: PetscCall(PetscMalloc1(npoints, &nnlist));
603: /* find nearest neighbour points in this cell */
604: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
605: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
606: for (q = 0; q < npoints; q++) {
607: PetscInt qn, nearest_neighbour = -1;
608: PetscReal sep, min_sep = PETSC_MAX_REAL;
610: coor_q = &xp[dim * q];
611: for (qn = 0; qn < npoints_e; qn++) {
612: coor_qn = &swarm_coor[dim * plist_e[qn]];
613: sep = 0.0;
614: for (d = 0; d < dim; d++) sep += (coor_q[d] - coor_qn[d]) * (coor_q[d] - coor_qn[d]);
615: if (sep < min_sep) {
616: nearest_neighbour = plist_e[qn];
617: min_sep = sep;
618: }
619: }
620: PetscCheck(nearest_neighbour != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell %" PetscInt_FMT " is empty - cannot initialize using nearest neighbours", e);
621: nnlist[q] = nearest_neighbour;
622: }
623: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
624: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
626: /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
627: for (q = 0; q < npoints; q++) PetscCall(DMSwarmCopyPoint(dm, nnlist[q], ncurr + q));
628: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
629: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
630: for (q = 0; q < npoints; q++) {
631: /* set the coordinates */
632: for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
633: /* set the cell index */
634: swarm_cellid[ncurr + q] = e;
635: }
636: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
637: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
639: PetscCall(DMSwarmSortRestorePointsPerCell(dm, e, &npoints_e, &plist_e));
640: PetscCall(PetscFree(nnlist));
641: } else {
642: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
643: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
644: for (q = 0; q < npoints; q++) {
645: /* set the coordinates */
646: for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
647: /* set the cell index */
648: swarm_cellid[ncurr + q] = e;
649: }
650: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
651: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
652: }
654: PetscCall(PetscFree(xp));
655: PetscCall(PetscFree(elcoor));
656: for (q = 0; q < npoints; q++) PetscCall(PetscFree(basis[q]));
657: PetscCall(PetscFree(basis));
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }
661: PetscErrorCode MaterialPoint_PopulateCell(DM dm_vp, DM dm_mpoint)
662: {
663: PetscInt _npe, _nel, e, nel;
664: const PetscInt *element;
665: DM dmc;
666: PetscQuadrature quadrature;
667: const PetscReal *xi;
668: PetscInt npoints_q, cnt;
670: PetscFunctionBeginUser;
671: PetscCall(DMDAGetElements(dm_vp, &_nel, &_npe, &element));
672: nel = _nel;
673: PetscCall(DMDARestoreElements(dm_vp, &_nel, &_npe, &element));
675: PetscCall(PetscDTGaussTensorQuadrature(2, 1, 4, -1.0, 1.0, &quadrature));
676: PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xi, NULL));
677: PetscCall(DMSwarmGetCellDM(dm_mpoint, &dmc));
679: PetscCall(DMSwarmSortGetAccess(dm_mpoint));
681: cnt = 0;
682: for (e = 0; e < nel; e++) {
683: PetscInt npoints_per_cell;
685: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm_mpoint, e, &npoints_per_cell));
687: if (npoints_per_cell < 12) {
688: PetscCall(DMSwarmPICInsertPointsCellwise(dm_mpoint, dm_vp, e, npoints_q, (PetscReal *)xi, PETSC_TRUE));
689: cnt++;
690: }
691: }
692: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &cnt, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
693: if (cnt > 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... ....pop cont: adjusted %" PetscInt_FMT " cells\n", cnt));
695: PetscCall(DMSwarmSortRestoreAccess(dm_mpoint));
696: PetscCall(PetscQuadratureDestroy(&quadrature));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: PetscErrorCode MaterialPoint_AdvectRK1(DM dm_vp, Vec vp, PetscReal dt, DM dm_mpoint)
701: {
702: DMSwarmCellDM celldm;
703: Vec vp_l, coor_l;
704: const PetscScalar *LA_vp;
705: PetscInt i, p, e, npoints, nel, npe, Nfc;
706: PetscInt *mpfield_cell;
707: PetscReal *mpfield_coor;
708: const PetscInt *element_list;
709: const PetscInt *element;
710: PetscScalar xi_p[NSD], Ni[NODES_PER_EL];
711: const PetscScalar *LA_coor;
712: PetscScalar dx[NSD];
713: const char **coordFields, *cellid;
715: PetscFunctionBeginUser;
716: PetscCall(DMGetCoordinatesLocal(dm_vp, &coor_l));
717: PetscCall(VecGetArrayRead(coor_l, &LA_coor));
719: PetscCall(DMGetLocalVector(dm_vp, &vp_l));
720: PetscCall(DMGlobalToLocalBegin(dm_vp, vp, INSERT_VALUES, vp_l));
721: PetscCall(DMGlobalToLocalEnd(dm_vp, vp, INSERT_VALUES, vp_l));
722: PetscCall(VecGetArrayRead(vp_l, &LA_vp));
724: PetscCall(DMDAGetElements(dm_vp, &nel, &npe, &element_list));
725: PetscCall(DMSwarmGetLocalSize(dm_mpoint, &npoints));
726: PetscCall(DMSwarmGetCellDMActive(dm_mpoint, &celldm));
727: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
728: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
729: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm_mpoint), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
730: PetscCall(DMSwarmGetField(dm_mpoint, coordFields[0], NULL, NULL, (void **)&mpfield_coor));
731: PetscCall(DMSwarmGetField(dm_mpoint, cellid, NULL, NULL, (void **)&mpfield_cell));
732: for (p = 0; p < npoints; p++) {
733: PetscReal *coor_p;
734: PetscScalar vel_n[NSD * NODES_PER_EL], vel_p[NSD];
735: const PetscScalar *x0;
736: const PetscScalar *x2;
738: e = mpfield_cell[p];
739: coor_p = &mpfield_coor[NSD * p];
740: element = &element_list[NODES_PER_EL * e];
742: /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
743: x0 = &LA_coor[NSD * element[0]];
744: x2 = &LA_coor[NSD * element[2]];
746: dx[0] = x2[0] - x0[0];
747: dx[1] = x2[1] - x0[1];
749: xi_p[0] = 2.0 * (coor_p[0] - x0[0]) / dx[0] - 1.0;
750: xi_p[1] = 2.0 * (coor_p[1] - x0[1]) / dx[1] - 1.0;
751: 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);
752: 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);
753: 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);
754: 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);
756: /* evaluate basis functions */
757: EvaluateBasis_Q1(xi_p, Ni);
759: /* get cell nodal velocities */
760: for (i = 0; i < NODES_PER_EL; i++) {
761: PetscInt nid;
763: nid = element[i];
764: vel_n[NSD * i + 0] = LA_vp[(NSD + 1) * nid + 0];
765: vel_n[NSD * i + 1] = LA_vp[(NSD + 1) * nid + 1];
766: }
768: /* interpolate velocity */
769: vel_p[0] = vel_p[1] = 0.0;
770: for (i = 0; i < NODES_PER_EL; i++) {
771: vel_p[0] += Ni[i] * vel_n[NSD * i + 0];
772: vel_p[1] += Ni[i] * vel_n[NSD * i + 1];
773: }
775: coor_p[0] += dt * PetscRealPart(vel_p[0]);
776: coor_p[1] += dt * PetscRealPart(vel_p[1]);
777: }
779: PetscCall(DMSwarmRestoreField(dm_mpoint, cellid, NULL, NULL, (void **)&mpfield_cell));
780: PetscCall(DMSwarmRestoreField(dm_mpoint, coordFields[0], NULL, NULL, (void **)&mpfield_coor));
781: PetscCall(DMDARestoreElements(dm_vp, &nel, &npe, &element_list));
782: PetscCall(VecRestoreArrayRead(vp_l, &LA_vp));
783: PetscCall(DMRestoreLocalVector(dm_vp, &vp_l));
784: PetscCall(VecRestoreArrayRead(coor_l, &LA_coor));
785: PetscFunctionReturn(PETSC_SUCCESS);
786: }
788: PetscErrorCode MaterialPoint_Interpolate(DM dm, Vec eta_v, Vec rho_v, DM dm_quadrature)
789: {
790: Vec eta_l, rho_l;
791: PetscScalar *_eta_l, *_rho_l;
792: PetscInt nqp, npe, nel;
793: PetscScalar qp_xi[GAUSS_POINTS][NSD];
794: PetscScalar qp_weight[GAUSS_POINTS];
795: PetscInt q, k, e;
796: PetscScalar Ni[GAUSS_POINTS][NODES_PER_EL];
797: const PetscInt *element_list;
798: PetscReal *q_eta, *q_rhs;
800: PetscFunctionBeginUser;
801: /* define quadrature rule */
802: CreateGaussQuadrature(&nqp, qp_xi, qp_weight);
803: for (q = 0; q < nqp; q++) EvaluateBasis_Q1(qp_xi[q], Ni[q]);
805: PetscCall(DMGetLocalVector(dm, &eta_l));
806: PetscCall(DMGetLocalVector(dm, &rho_l));
808: PetscCall(DMGlobalToLocalBegin(dm, eta_v, INSERT_VALUES, eta_l));
809: PetscCall(DMGlobalToLocalEnd(dm, eta_v, INSERT_VALUES, eta_l));
810: PetscCall(DMGlobalToLocalBegin(dm, rho_v, INSERT_VALUES, rho_l));
811: PetscCall(DMGlobalToLocalEnd(dm, rho_v, INSERT_VALUES, rho_l));
813: PetscCall(VecGetArray(eta_l, &_eta_l));
814: PetscCall(VecGetArray(rho_l, &_rho_l));
816: PetscCall(DMSwarmGetField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
817: PetscCall(DMSwarmGetField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
819: PetscCall(DMDAGetElements(dm, &nel, &npe, &element_list));
820: for (e = 0; e < nel; e++) {
821: PetscScalar eta_field_e[NODES_PER_EL];
822: PetscScalar rho_field_e[NODES_PER_EL];
823: const PetscInt *element = &element_list[4 * e];
825: for (k = 0; k < NODES_PER_EL; k++) {
826: eta_field_e[k] = _eta_l[element[k]];
827: rho_field_e[k] = _rho_l[element[k]];
828: }
830: for (q = 0; q < nqp; q++) {
831: PetscScalar eta_q, rho_q;
833: eta_q = rho_q = 0.0;
834: for (k = 0; k < NODES_PER_EL; k++) {
835: eta_q += Ni[q][k] * eta_field_e[k];
836: rho_q += Ni[q][k] * rho_field_e[k];
837: }
839: q_eta[nqp * e + q] = PetscRealPart(eta_q);
840: q_rhs[nqp * e + q] = PetscRealPart(rho_q);
841: }
842: }
843: PetscCall(DMDARestoreElements(dm, &nel, &npe, &element_list));
845: PetscCall(DMSwarmRestoreField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
846: PetscCall(DMSwarmRestoreField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
848: PetscCall(VecRestoreArray(rho_l, &_rho_l));
849: PetscCall(VecRestoreArray(eta_l, &_eta_l));
850: PetscCall(DMRestoreLocalVector(dm, &rho_l));
851: PetscCall(DMRestoreLocalVector(dm, &eta_l));
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: static PetscErrorCode SolveTimeDepStokes(PetscInt mx, PetscInt my)
856: {
857: DM dm_stokes, dm_coeff;
858: PetscInt u_dof, p_dof, dof, stencil_width;
859: Mat A, B;
860: PetscInt nel_local;
861: Vec eta_v, rho_v;
862: Vec f, X;
863: KSP ksp;
864: PC pc;
865: char filename[PETSC_MAX_PATH_LEN];
866: DM dms_quadrature, dms_mpoint;
867: PetscInt nel, npe, npoints;
868: const PetscInt *element_list;
869: PetscInt tk, nt, dump_freq;
870: PetscReal dt, dt_max = 0.0;
871: PetscReal vx[2], vy[2], max_v = 0.0, max_v_step, dh;
872: const char *fieldnames[] = {"eta", "rho"};
873: Vec pfields[2];
874: PetscInt ppcell = 1;
875: PetscReal time, delta_eta = 1.0;
876: PetscBool randomize_coords = PETSC_FALSE;
877: PetscReal randomize_fac = 0.25;
878: PetscBool no_view = PETSC_FALSE;
879: PetscBool isbddc;
881: PetscFunctionBeginUser;
882: /*
883: Generate the DMDA for the velocity and pressure spaces.
884: We use Q1 elements for both fields.
885: The Q1 FE basis on a regular mesh has a 9-point stencil (DMDA_STENCIL_BOX)
886: The number of nodes in each direction is mx+1, my+1
887: */
888: u_dof = U_DOFS; /* Vx, Vy - velocities */
889: p_dof = P_DOFS; /* p - pressure */
890: dof = u_dof + p_dof;
891: stencil_width = 1;
892: 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));
893: PetscCall(DMDASetElementType(dm_stokes, DMDA_ELEMENT_Q1));
894: PetscCall(DMSetMatType(dm_stokes, MATAIJ));
895: PetscCall(DMSetFromOptions(dm_stokes));
896: PetscCall(DMSetUp(dm_stokes));
897: PetscCall(DMDASetFieldName(dm_stokes, 0, "ux"));
898: PetscCall(DMDASetFieldName(dm_stokes, 1, "uy"));
899: PetscCall(DMDASetFieldName(dm_stokes, 2, "p"));
901: /* unit box [0,0.9142] x [0,1] */
902: PetscCall(DMDASetUniformCoordinates(dm_stokes, 0.0, 0.9142, 0.0, 1.0, 0., 0.));
903: dh = 1.0 / (PetscReal)mx;
905: /* Get local number of elements */
906: {
907: PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));
909: nel_local = nel;
911: PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
912: }
914: /* Create DMDA for representing scalar fields */
915: PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 1, &dm_coeff));
917: /* Create the swarm for storing quadrature point values */
918: PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_quadrature));
919: PetscCall(DMSetType(dms_quadrature, DMSWARM));
920: PetscCall(DMSetDimension(dms_quadrature, 2));
921: PetscCall(PetscObjectSetName((PetscObject)dms_quadrature, "Quadrature Swarm"));
923: /* Register fields for viscosity and density on the quadrature points */
924: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "eta_q", 1, PETSC_REAL));
925: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "rho_q", 1, PETSC_REAL));
926: PetscCall(DMSwarmFinalizeFieldRegister(dms_quadrature));
927: PetscCall(DMSwarmSetLocalSizes(dms_quadrature, nel_local * GAUSS_POINTS, 0));
929: /* Create the material point swarm */
930: PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_mpoint));
931: PetscCall(DMSetType(dms_mpoint, DMSWARM));
932: PetscCall(DMSetDimension(dms_mpoint, 2));
933: PetscCall(PetscObjectSetName((PetscObject)dms_mpoint, "Material Point Swarm"));
935: /* Configure the material point swarm to be of type Particle-In-Cell */
936: PetscCall(DMSwarmSetType(dms_mpoint, DMSWARM_PIC));
938: /*
939: Specify the DM to use for point location and projections
940: within the context of a PIC scheme
941: */
942: PetscCall(DMSwarmSetCellDM(dms_mpoint, dm_coeff));
944: /* Register fields for viscosity and density */
945: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "eta", 1, PETSC_REAL));
946: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "rho", 1, PETSC_REAL));
947: PetscCall(DMSwarmFinalizeFieldRegister(dms_mpoint));
949: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL));
950: PetscCall(DMSwarmSetLocalSizes(dms_mpoint, nel_local * ppcell, 100));
952: /*
953: Layout the material points in space using the cell DM.
954: Particle coordinates are defined by cell wise using different methods.
955: - DMSWARMPIC_LAYOUT_GAUSS defines particles coordinates at the positions
956: corresponding to a Gauss quadrature rule with
957: ppcell points in each direction.
958: - DMSWARMPIC_LAYOUT_REGULAR defines particle coordinates at the centoid of
959: ppcell x ppcell quadralaterals defined within the
960: reference element.
961: - DMSWARMPIC_LAYOUT_SUBDIVISION defines particles coordinates at the centroid
962: of each quadralateral obtained by sub-dividing
963: the reference element cell ppcell times.
964: */
965: PetscCall(DMSwarmInsertPointsUsingCellDM(dms_mpoint, DMSWARMPIC_LAYOUT_SUBDIVISION, ppcell));
967: /*
968: Defne a high resolution layer of material points across the material interface
969: */
970: {
971: PetscInt npoints_dir_x[2];
972: PetscReal min[2], max[2];
974: npoints_dir_x[0] = (PetscInt)(0.9142 / (0.05 * dh));
975: npoints_dir_x[1] = (PetscInt)((0.25 - 0.15) / (0.05 * dh));
976: min[0] = 0.0;
977: max[0] = 0.9142;
978: min[1] = 0.05;
979: max[1] = 0.35;
980: PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
981: }
983: /*
984: Define a high resolution layer of material points near the surface of the domain
985: to deal with weakly compressible Q1-Q1 elements. These elements "self compact"
986: when applied to buoyancy driven flow. The error in div(u) is O(h).
987: */
988: {
989: PetscInt npoints_dir_x[2];
990: PetscReal min[2], max[2];
992: npoints_dir_x[0] = (PetscInt)(0.9142 / (0.25 * dh));
993: npoints_dir_x[1] = (PetscInt)(3.0 * dh / (0.25 * dh));
994: min[0] = 0.0;
995: max[0] = 0.9142;
996: min[1] = 1.0 - 3.0 * dh;
997: max[1] = 1.0 - 0.0001;
998: PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
999: }
1001: PetscCall(DMView(dms_mpoint, PETSC_VIEWER_STDOUT_WORLD));
1003: /* Define initial material properties on each particle in the material point swarm */
1004: PetscCall(PetscOptionsGetReal(NULL, NULL, "-delta_eta", &delta_eta, NULL));
1005: PetscCall(PetscOptionsGetBool(NULL, NULL, "-randomize_coords", &randomize_coords, NULL));
1006: PetscCall(PetscOptionsGetReal(NULL, NULL, "-randomize_fac", &randomize_fac, NULL));
1007: PetscCheck(randomize_fac <= 1.0, PETSC_COMM_WORLD, PETSC_ERR_USER, "The value of -randomize_fac should be <= 1.0");
1008: {
1009: PetscReal *array_x, *array_e, *array_r;
1010: PetscRandom r;
1011: PetscMPIInt rank;
1013: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1015: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
1016: PetscCall(PetscRandomSetInterval(r, -randomize_fac * dh, randomize_fac * dh));
1017: PetscCall(PetscRandomSetSeed(r, (unsigned long)rank));
1018: PetscCall(PetscRandomSeed(r));
1020: PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));
1022: /*
1023: Fetch the registered data from the material point DMSwarm.
1024: The fields "eta" and "rho" were registered by this example.
1025: The field identified by the variable DMSwarmPICField_coor
1026: was registered by the DMSwarm implementation when the function
1027: DMSwarmSetType(dms_mpoint,DMSWARM_PIC)
1028: was called. The returned array defines the coordinates of each
1029: material point in the point swarm.
1030: */
1031: PetscCall(DMSwarmGetField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));
1032: PetscCall(DMSwarmGetField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
1033: PetscCall(DMSwarmGetField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));
1035: PetscCall(DMSwarmGetLocalSize(dms_mpoint, &npoints));
1036: for (PetscInt p = 0; p < npoints; p++) {
1037: PetscReal x_p[2], rr[2];
1039: if (randomize_coords) {
1040: PetscCall(PetscRandomGetValueReal(r, &rr[0]));
1041: PetscCall(PetscRandomGetValueReal(r, &rr[1]));
1042: array_x[2 * p + 0] += rr[0];
1043: array_x[2 * p + 1] += rr[1];
1044: }
1046: /* Get the coordinates of point, p */
1047: x_p[0] = array_x[2 * p + 0];
1048: x_p[1] = array_x[2 * p + 1];
1050: if (x_p[1] < (0.2 + 0.02 * PetscCosReal(PETSC_PI * x_p[0] / 0.9142))) {
1051: /* Material properties below the interface */
1052: array_e[p] = 1.0 * (1.0 / delta_eta);
1053: array_r[p] = 0.0;
1054: } else {
1055: /* Material properties above the interface */
1056: array_e[p] = 1.0;
1057: array_r[p] = 1.0;
1058: }
1059: }
1061: /*
1062: Restore the fetched data fields from the material point DMSwarm.
1063: Calling the Restore function invalidates the points array_r, array_e, array_x
1064: by setting them to NULL.
1065: */
1066: PetscCall(DMSwarmRestoreField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));
1067: PetscCall(DMSwarmRestoreField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
1068: PetscCall(DMSwarmRestoreField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));
1070: PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
1071: PetscCall(PetscRandomDestroy(&r));
1072: }
1074: /*
1075: If the particle coordinates where randomly shifted, they may have crossed into another
1076: element, or into another sub-domain. To account for this we call the Migrate function.
1077: */
1078: if (randomize_coords) PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));
1080: PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL));
1081: if (!no_view) PetscCall(DMSwarmViewXDMF(dms_mpoint, "ic_coeff_dms.xmf"));
1083: /* project the swarm properties */
1084: PetscCall(DMCreateGlobalVector(dm_coeff, &pfields[0]));
1085: PetscCall(DMCreateGlobalVector(dm_coeff, &pfields[1]));
1086: PetscCall(DMSwarmProjectFields(dms_mpoint, NULL, 2, fieldnames, pfields, SCATTER_FORWARD));
1087: eta_v = pfields[0];
1088: rho_v = pfields[1];
1089: PetscCall(PetscObjectSetName((PetscObject)eta_v, "eta"));
1090: PetscCall(PetscObjectSetName((PetscObject)rho_v, "rho"));
1091: PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));
1093: /* view projected coefficients eta and rho */
1094: if (!no_view) {
1095: PetscViewer viewer;
1097: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1098: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1099: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1100: PetscCall(PetscViewerFileSetName(viewer, "ic_coeff_dmda.vts"));
1101: PetscCall(VecView(eta_v, viewer));
1102: PetscCall(VecView(rho_v, viewer));
1103: PetscCall(PetscViewerDestroy(&viewer));
1104: }
1106: PetscCall(DMCreateMatrix(dm_stokes, &A));
1107: PetscCall(DMCreateMatrix(dm_stokes, &B));
1108: PetscCall(DMCreateGlobalVector(dm_stokes, &f));
1109: PetscCall(DMCreateGlobalVector(dm_stokes, &X));
1111: PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
1112: PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
1113: PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));
1115: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
1116: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));
1118: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
1119: PetscCall(KSPSetOptionsPrefix(ksp, "stokes_"));
1120: PetscCall(KSPSetDM(ksp, dm_stokes));
1121: PetscCall(KSPSetDMActive(ksp, KSP_DMACTIVE_ALL, PETSC_FALSE));
1122: PetscCall(KSPSetOperators(ksp, A, B));
1123: PetscCall(KSPSetFromOptions(ksp));
1124: PetscCall(KSPGetPC(ksp, &pc));
1125: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &isbddc));
1126: if (isbddc) PetscCall(KSPSetOperators(ksp, A, A));
1128: /* Define u-v-p indices for fieldsplit */
1129: {
1130: PC pc;
1131: const PetscInt ufields[] = {0, 1}, pfields[1] = {2};
1133: PetscCall(KSPGetPC(ksp, &pc));
1134: PetscCall(PCFieldSplitSetBlockSize(pc, 3));
1135: PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
1136: PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
1137: }
1139: /* If using a fieldsplit preconditioner, attach a DMDA to the velocity split so that geometric multigrid can be used */
1140: {
1141: PC pc, pc_u;
1142: KSP *sub_ksp, ksp_u;
1143: PetscInt nsplits;
1144: DM dm_u;
1145: PetscBool is_pcfs;
1147: PetscCall(KSPGetPC(ksp, &pc));
1149: is_pcfs = PETSC_FALSE;
1150: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs));
1152: if (is_pcfs) {
1153: PetscCall(KSPSetUp(ksp));
1154: PetscCall(KSPGetPC(ksp, &pc));
1155: PetscCall(PCFieldSplitGetSubKSP(pc, &nsplits, &sub_ksp));
1156: ksp_u = sub_ksp[0];
1157: PetscCall(PetscFree(sub_ksp));
1159: if (nsplits == 2) {
1160: PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 2, &dm_u));
1162: PetscCall(KSPSetDM(ksp_u, dm_u));
1163: PetscCall(KSPSetDMActive(ksp_u, KSP_DMACTIVE_ALL, PETSC_FALSE));
1164: PetscCall(DMDestroy(&dm_u));
1166: /* enforce galerkin coarse grids be used */
1167: PetscCall(KSPGetPC(ksp_u, &pc_u));
1168: PetscCall(PCMGSetGalerkin(pc_u, PC_MG_GALERKIN_PMAT));
1169: }
1170: }
1171: }
1173: dump_freq = 10;
1174: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dump_freq", &dump_freq, NULL));
1175: nt = 10;
1176: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nt", &nt, NULL));
1177: time = 0.0;
1178: for (tk = 1; tk <= nt; tk++) {
1179: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... assemble\n"));
1180: PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
1181: PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
1182: PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));
1184: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... bc imposition\n"));
1185: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
1186: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));
1188: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... solve\n"));
1189: PetscCall(KSPSetOperators(ksp, A, isbddc ? A : B));
1190: PetscCall(KSPSolve(ksp, f, X));
1192: PetscCall(VecStrideMax(X, 0, NULL, &vx[1]));
1193: PetscCall(VecStrideMax(X, 1, NULL, &vy[1]));
1194: PetscCall(VecStrideMin(X, 0, NULL, &vx[0]));
1195: PetscCall(VecStrideMin(X, 1, NULL, &vy[0]));
1197: max_v_step = PetscMax(vx[0], vx[1]);
1198: max_v_step = PetscMax(max_v_step, vy[0]);
1199: max_v_step = PetscMax(max_v_step, vy[1]);
1200: max_v = PetscMax(max_v, max_v_step);
1202: dt_max = 2.0;
1203: dt = 0.5 * (dh / max_v_step);
1204: 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));
1205: dt = PetscMin(dt_max, dt);
1207: /* advect */
1208: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... advect\n"));
1209: PetscCall(MaterialPoint_AdvectRK1(dm_stokes, X, dt, dms_mpoint));
1211: /* migrate */
1212: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... migrate\n"));
1213: PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));
1215: /* update cell population */
1216: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... populate cells\n"));
1217: PetscCall(MaterialPoint_PopulateCell(dm_stokes, dms_mpoint));
1219: /* update coefficients on quadrature points */
1220: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... project\n"));
1221: PetscCall(DMSwarmProjectFields(dms_mpoint, NULL, 2, fieldnames, pfields, SCATTER_FORWARD));
1222: eta_v = pfields[0];
1223: rho_v = pfields[1];
1224: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... interp\n"));
1225: PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));
1227: if (tk % dump_freq == 0) {
1228: PetscViewer viewer;
1230: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... write XDMF, VTS\n"));
1231: PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_coeff_dms.xmf", tk));
1232: PetscCall(DMSwarmViewXDMF(dms_mpoint, filename));
1234: PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_vp_dm.vts", tk));
1235: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1236: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1237: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1238: PetscCall(PetscViewerFileSetName(viewer, filename));
1239: PetscCall(VecView(X, viewer));
1240: PetscCall(PetscViewerDestroy(&viewer));
1241: }
1242: time += dt;
1243: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step %" PetscInt_FMT " : time %1.2e \n", tk, (double)time));
1244: }
1246: PetscCall(KSPDestroy(&ksp));
1247: PetscCall(VecDestroy(&X));
1248: PetscCall(VecDestroy(&f));
1249: PetscCall(MatDestroy(&A));
1250: PetscCall(MatDestroy(&B));
1251: PetscCall(VecDestroy(&eta_v));
1252: PetscCall(VecDestroy(&rho_v));
1254: PetscCall(DMDestroy(&dms_mpoint));
1255: PetscCall(DMDestroy(&dms_quadrature));
1256: PetscCall(DMDestroy(&dm_coeff));
1257: PetscCall(DMDestroy(&dm_stokes));
1258: PetscFunctionReturn(PETSC_SUCCESS);
1259: }
1261: /*
1262: <sequential run>
1263: ./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
1264: */
1265: int main(int argc, char **args)
1266: {
1267: PetscInt mx, my;
1268: PetscBool set = PETSC_FALSE;
1270: PetscFunctionBeginUser;
1271: PetscCall(PetscInitialize(&argc, &args, NULL, help));
1272: mx = my = 10;
1273: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1274: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1275: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mxy", &mx, &set));
1276: if (set) my = mx;
1277: PetscCall(SolveTimeDepStokes(mx, my));
1278: PetscCall(PetscFinalize());
1279: return 0;
1280: }
1282: /* -------------------------- helpers for boundary conditions -------------------------------- */
1283: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
1284: {
1285: DM cda;
1286: Vec coords;
1287: PetscInt si, sj, nx, ny, i, j;
1288: PetscInt M, N;
1289: DMDACoor2d **_coords;
1290: const PetscInt *g_idx;
1291: PetscInt *bc_global_ids;
1292: PetscScalar *bc_vals;
1293: PetscInt nbcs;
1294: PetscInt n_dofs;
1295: ISLocalToGlobalMapping ltogm;
1297: PetscFunctionBeginUser;
1298: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1299: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1301: PetscCall(DMGetCoordinateDM(da, &cda));
1302: PetscCall(DMGetCoordinatesLocal(da, &coords));
1303: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1304: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1305: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1307: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1308: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1310: /* init the entries to -1 so VecSetValues will ignore them */
1311: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1313: i = nx - 1;
1314: for (j = 0; j < ny; j++) {
1315: PetscInt local_id;
1317: local_id = i + j * nx;
1319: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1321: bc_vals[j] = 0.0;
1322: }
1323: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1324: nbcs = 0;
1325: if ((si + nx) == (M)) nbcs = ny;
1327: if (b) {
1328: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1329: PetscCall(VecAssemblyBegin(b));
1330: PetscCall(VecAssemblyEnd(b));
1331: }
1332: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1334: PetscCall(PetscFree(bc_vals));
1335: PetscCall(PetscFree(bc_global_ids));
1337: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1338: PetscFunctionReturn(PETSC_SUCCESS);
1339: }
1341: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
1342: {
1343: DM cda;
1344: Vec coords;
1345: PetscInt si, sj, nx, ny, i, j;
1346: PetscInt M, N;
1347: DMDACoor2d **_coords;
1348: const PetscInt *g_idx;
1349: PetscInt *bc_global_ids;
1350: PetscScalar *bc_vals;
1351: PetscInt nbcs;
1352: PetscInt n_dofs;
1353: ISLocalToGlobalMapping ltogm;
1355: PetscFunctionBeginUser;
1356: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1357: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1359: PetscCall(DMGetCoordinateDM(da, &cda));
1360: PetscCall(DMGetCoordinatesLocal(da, &coords));
1361: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1362: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1363: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1365: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1366: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1368: /* init the entries to -1 so VecSetValues will ignore them */
1369: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1371: i = 0;
1372: for (j = 0; j < ny; j++) {
1373: PetscInt local_id;
1375: local_id = i + j * nx;
1377: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1379: bc_vals[j] = 0.0;
1380: }
1381: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1382: nbcs = 0;
1383: if (si == 0) nbcs = ny;
1385: if (b) {
1386: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1387: PetscCall(VecAssemblyBegin(b));
1388: PetscCall(VecAssemblyEnd(b));
1389: }
1391: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1393: PetscCall(PetscFree(bc_vals));
1394: PetscCall(PetscFree(bc_global_ids));
1396: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1397: PetscFunctionReturn(PETSC_SUCCESS);
1398: }
1400: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
1401: {
1402: DM cda;
1403: Vec coords;
1404: PetscInt si, sj, nx, ny, i, j;
1405: PetscInt M, N;
1406: DMDACoor2d **_coords;
1407: const PetscInt *g_idx;
1408: PetscInt *bc_global_ids;
1409: PetscScalar *bc_vals;
1410: PetscInt nbcs;
1411: PetscInt n_dofs;
1412: ISLocalToGlobalMapping ltogm;
1414: PetscFunctionBeginUser;
1415: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1416: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1418: PetscCall(DMGetCoordinateDM(da, &cda));
1419: PetscCall(DMGetCoordinatesLocal(da, &coords));
1420: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1421: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1422: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1424: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1425: PetscCall(PetscMalloc1(nx, &bc_vals));
1427: /* init the entries to -1 so VecSetValues will ignore them */
1428: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1430: j = ny - 1;
1431: for (i = 0; i < nx; i++) {
1432: PetscInt local_id;
1434: local_id = i + j * nx;
1436: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1438: bc_vals[i] = 0.0;
1439: }
1440: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1441: nbcs = 0;
1442: if ((sj + ny) == (N)) nbcs = nx;
1444: if (b) {
1445: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1446: PetscCall(VecAssemblyBegin(b));
1447: PetscCall(VecAssemblyEnd(b));
1448: }
1449: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, NULL, NULL));
1451: PetscCall(PetscFree(bc_vals));
1452: PetscCall(PetscFree(bc_global_ids));
1454: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1455: PetscFunctionReturn(PETSC_SUCCESS);
1456: }
1458: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
1459: {
1460: DM cda;
1461: Vec coords;
1462: PetscInt si, sj, nx, ny, i, j;
1463: PetscInt M, N;
1464: DMDACoor2d **_coords;
1465: const PetscInt *g_idx;
1466: PetscInt *bc_global_ids;
1467: PetscScalar *bc_vals;
1468: PetscInt nbcs;
1469: PetscInt n_dofs;
1470: ISLocalToGlobalMapping ltogm;
1472: PetscFunctionBeginUser;
1473: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1474: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1476: PetscCall(DMGetCoordinateDM(da, &cda));
1477: PetscCall(DMGetCoordinatesLocal(da, &coords));
1478: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1479: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1480: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1482: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1483: PetscCall(PetscMalloc1(nx, &bc_vals));
1485: /* init the entries to -1 so VecSetValues will ignore them */
1486: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1488: j = 0;
1489: for (i = 0; i < nx; i++) {
1490: PetscInt local_id;
1492: local_id = i + j * nx;
1494: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1496: bc_vals[i] = 0.0;
1497: }
1498: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1499: nbcs = 0;
1500: if (sj == 0) nbcs = nx;
1502: if (b) {
1503: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1504: PetscCall(VecAssemblyBegin(b));
1505: PetscCall(VecAssemblyEnd(b));
1506: }
1507: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1509: PetscCall(PetscFree(bc_vals));
1510: PetscCall(PetscFree(bc_global_ids));
1512: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1513: PetscFunctionReturn(PETSC_SUCCESS);
1514: }
1516: /*
1517: Impose free slip boundary conditions on the left/right faces: u_i n_i = 0, tau_{ij} t_j = 0
1518: Impose no slip boundray conditions on the top/bottom faces: u_i n_i = 0, u_i t_i = 0
1519: */
1520: static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes, Mat A, Vec f)
1521: {
1522: PetscFunctionBeginUser;
1523: PetscCall(BCApplyZero_NORTH(dm_stokes, 0, A, f));
1524: PetscCall(BCApplyZero_NORTH(dm_stokes, 1, A, f));
1525: PetscCall(BCApplyZero_EAST(dm_stokes, 0, A, f));
1526: PetscCall(BCApplyZero_SOUTH(dm_stokes, 0, A, f));
1527: PetscCall(BCApplyZero_SOUTH(dm_stokes, 1, A, f));
1528: PetscCall(BCApplyZero_WEST(dm_stokes, 0, A, f));
1529: PetscFunctionReturn(PETSC_SUCCESS);
1530: }
1532: /*TEST
1534: test:
1535: suffix: 1
1536: args: -no_view
1537: requires: !complex double
1538: filter: grep -v atomic
1539: filter_output: grep -v atomic
1540: test:
1541: suffix: 1_matis
1542: requires: !complex double
1543: args: -no_view -dm_mat_type is
1544: filter: grep -v atomic
1545: filter_output: grep -v atomic
1546: testset:
1547: nsize: 4
1548: requires: !complex double
1549: 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
1550: filter: grep -v atomic
1551: filter_output: grep -v atomic
1552: test:
1553: suffix: fetidp
1554: args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
1555: test:
1556: suffix: fetidp_lumped
1557: 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
1558: test:
1559: suffix: fetidp_saddlepoint
1560: 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
1561: test:
1562: suffix: fetidp_saddlepoint_lumped
1563: 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
1564: TEST*/