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>
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: PetscFunctionBeginUser;
122: for (PetscInt i = 0; i < NODES_PER_EL; i++) {
123: /* velocity */
124: s_u[NSD * i + 0] = 3 * element[i];
125: s_u[NSD * i + 1] = 3 * element[i] + 1;
126: /* pressure */
127: s_p[i] = 3 * element[i] + 2;
128: }
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: 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)
133: {
134: PetscInt ij, r, c, nc;
136: nc = u_NPE * u_dof;
137: r = w_dof * wi + wd;
138: c = u_dof * ui + ud;
139: ij = r * nc + c;
140: return ij;
141: }
143: static void BForm_DivT(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
144: {
145: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
146: PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
147: PetscScalar J_p, tildeD[3];
148: PetscScalar B[3][U_DOFS * NODES_PER_EL];
149: PetscInt p, i, j, k, ngp;
151: /* define quadrature rule */
152: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
154: /* evaluate bilinear form */
155: for (p = 0; p < ngp; p++) {
156: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
157: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
159: for (i = 0; i < NODES_PER_EL; i++) {
160: PetscScalar d_dx_i = GNx_p[0][i];
161: PetscScalar d_dy_i = GNx_p[1][i];
163: B[0][2 * i] = d_dx_i;
164: B[0][2 * i + 1] = 0.0;
165: B[1][2 * i] = 0.0;
166: B[1][2 * i + 1] = d_dy_i;
167: B[2][2 * i] = d_dy_i;
168: B[2][2 * i + 1] = d_dx_i;
169: }
171: tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
172: tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
173: tildeD[2] = gp_weight[p] * J_p * eta[p];
175: /* form Bt tildeD B */
176: /*
177: Ke_ij = Bt_ik . D_kl . B_lj
178: = B_ki . D_kl . B_lj
179: = B_ki . D_kk . B_kj
180: */
181: for (i = 0; i < 8; i++) {
182: for (j = 0; j < 8; j++) {
183: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
184: Ke[i + 8 * j] += B[k][i] * tildeD[k] * B[k][j];
185: }
186: }
187: }
188: }
189: }
191: static void BForm_Grad(PetscScalar Ke[], PetscScalar coords[])
192: {
193: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
194: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
195: PetscScalar J_p, fac;
196: PetscInt p, i, j, di, ngp;
198: /* define quadrature rule */
199: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
201: /* evaluate bilinear form */
202: for (p = 0; p < ngp; p++) {
203: EvaluateBasis_Q1(gp_xi[p], Ni_p);
204: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
205: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
206: fac = gp_weight[p] * J_p;
208: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
209: for (di = 0; di < NSD; di++) { /* u dofs */
210: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
211: PetscInt IJ;
212: IJ = map_wIwDI_uJuDJ(i, di, NODES_PER_EL, 2, j, 0, NODES_PER_EL, 1);
214: Ke[IJ] -= GNx_p[di][i] * Ni_p[j] * fac;
215: }
216: }
217: }
218: }
219: }
221: static void BForm_Div(PetscScalar De[], PetscScalar coords[])
222: {
223: PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
224: PetscInt i, j, nr_g, nc_g;
226: PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));
227: BForm_Grad(Ge, coords);
229: nr_g = U_DOFS * NODES_PER_EL;
230: nc_g = P_DOFS * NODES_PER_EL;
232: for (i = 0; i < nr_g; i++) {
233: for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
234: }
235: }
237: static void BForm_Stabilisation(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
238: {
239: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
240: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
241: PetscScalar J_p, fac, eta_avg;
242: PetscInt p, i, j, ngp;
244: /* define quadrature rule */
245: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
247: /* evaluate bilinear form */
248: for (p = 0; p < ngp; p++) {
249: EvaluateBasis_Q1(gp_xi[p], Ni_p);
250: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
251: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
252: fac = gp_weight[p] * J_p;
254: for (i = 0; i < NODES_PER_EL; i++) {
255: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * (Ni_p[i] * Ni_p[j] - 0.0625);
256: }
257: }
259: /* scale */
260: eta_avg = 0.0;
261: for (p = 0; p < ngp; p++) eta_avg += eta[p];
262: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
263: fac = 1.0 / eta_avg;
264: for (i = 0; i < NODES_PER_EL; i++) {
265: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
266: }
267: }
269: static void BForm_ScaledMassMatrix(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
270: {
271: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
272: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
273: PetscScalar J_p, fac, eta_avg;
274: PetscInt p, i, j, ngp;
276: /* define quadrature rule */
277: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
279: /* evaluate bilinear form */
280: for (p = 0; p < ngp; p++) {
281: EvaluateBasis_Q1(gp_xi[p], Ni_p);
282: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
283: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
284: fac = gp_weight[p] * J_p;
286: for (i = 0; i < NODES_PER_EL; i++) {
287: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * Ni_p[i] * Ni_p[j];
288: }
289: }
291: /* scale */
292: eta_avg = 0.0;
293: for (p = 0; p < ngp; p++) eta_avg += eta[p];
294: eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
295: fac = 1.0 / eta_avg;
296: for (i = 0; i < NODES_PER_EL; i++) {
297: for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] *= fac;
298: }
299: }
301: static void LForm_MomentumRHS(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
302: {
303: PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
304: PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
305: PetscScalar J_p, fac;
306: PetscInt p, i, ngp;
308: /* define quadrature rule */
309: CreateGaussQuadrature(&ngp, gp_xi, gp_weight);
311: /* evaluate linear form */
312: for (p = 0; p < ngp; p++) {
313: EvaluateBasis_Q1(gp_xi[p], Ni_p);
314: EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
315: EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
316: fac = gp_weight[p] * J_p;
318: for (i = 0; i < NODES_PER_EL; i++) {
319: Fe[NSD * i] = 0.0;
320: Fe[NSD * i + 1] -= fac * Ni_p[i] * fy[p];
321: }
322: }
323: }
325: static PetscErrorCode GetElementCoords(const PetscScalar _coords[], const PetscInt e2n[], PetscScalar el_coords[])
326: {
327: PetscFunctionBeginUser;
328: /* get coords for the element */
329: for (PetscInt i = 0; i < 4; i++) {
330: for (PetscInt d = 0; d < NSD; d++) el_coords[NSD * i + d] = _coords[NSD * e2n[i] + d];
331: }
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: static PetscErrorCode AssembleStokes_A(Mat A, DM stokes_da, DM quadrature)
336: {
337: DM cda;
338: Vec coords;
339: const PetscScalar *_coords;
340: PetscInt u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
341: PetscInt p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
342: PetscInt nel, npe, eidx;
343: const PetscInt *element_list;
344: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
345: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
346: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
347: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
348: PetscScalar el_coords[NODES_PER_EL * NSD];
349: PetscScalar *q_eta, *prop_eta;
351: PetscFunctionBeginUser;
352: PetscCall(MatZeroEntries(A));
353: /* setup for coords */
354: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
355: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
356: PetscCall(VecGetArrayRead(coords, &_coords));
358: /* setup for coefficients */
359: PetscCall(DMSwarmGetField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
361: PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
362: for (eidx = 0; eidx < nel; eidx++) {
363: const PetscInt *element = &element_list[npe * eidx];
365: /* get coords for the element */
366: PetscCall(GetElementCoords(_coords, element, el_coords));
368: /* get coefficients for the element */
369: prop_eta = &q_eta[GAUSS_POINTS * eidx];
371: /* initialise element stiffness matrix */
372: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
373: PetscCall(PetscMemzero(Ge, sizeof(Ge)));
374: PetscCall(PetscMemzero(De, sizeof(De)));
375: PetscCall(PetscMemzero(Ce, sizeof(Ce)));
377: /* form element stiffness matrix */
378: BForm_DivT(Ae, el_coords, prop_eta);
379: BForm_Grad(Ge, el_coords);
380: BForm_Div(De, el_coords);
381: BForm_Stabilisation(Ce, el_coords, prop_eta);
383: /* insert element matrix into global matrix */
384: PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
385: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
386: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
387: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
388: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
389: }
390: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
391: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
393: PetscCall(DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
394: PetscCall(VecRestoreArrayRead(coords, &_coords));
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: static PetscErrorCode AssembleStokes_PC(Mat A, DM stokes_da, DM quadrature)
399: {
400: DM cda;
401: Vec coords;
402: const PetscScalar *_coords;
403: PetscInt u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
404: PetscInt p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
405: PetscInt nel, npe, eidx;
406: const PetscInt *element_list;
407: PetscScalar Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
408: PetscScalar Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
409: PetscScalar De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
410: PetscScalar Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
411: PetscScalar el_coords[NODES_PER_EL * NSD];
412: PetscScalar *q_eta, *prop_eta;
414: PetscFunctionBeginUser;
415: PetscCall(MatZeroEntries(A));
416: /* setup for coords */
417: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
418: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
419: PetscCall(VecGetArrayRead(coords, &_coords));
421: /* setup for coefficients */
422: PetscCall(DMSwarmGetField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
424: PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
425: for (eidx = 0; eidx < nel; eidx++) {
426: const PetscInt *element = &element_list[npe * eidx];
428: /* get coords for the element */
429: PetscCall(GetElementCoords(_coords, element, el_coords));
431: /* get coefficients for the element */
432: prop_eta = &q_eta[GAUSS_POINTS * eidx];
434: /* initialise element stiffness matrix */
435: PetscCall(PetscMemzero(Ae, sizeof(Ae)));
436: PetscCall(PetscMemzero(Ge, sizeof(Ge)));
437: PetscCall(PetscMemzero(De, sizeof(De)));
438: PetscCall(PetscMemzero(Ce, sizeof(Ce)));
440: /* form element stiffness matrix */
441: BForm_DivT(Ae, el_coords, prop_eta);
442: BForm_Grad(Ge, el_coords);
443: BForm_ScaledMassMatrix(Ce, el_coords, prop_eta);
445: /* insert element matrix into global matrix */
446: PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
447: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
448: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
449: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
450: PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
451: }
452: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
453: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
455: PetscCall(DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
456: PetscCall(VecRestoreArrayRead(coords, &_coords));
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: static PetscErrorCode AssembleStokes_RHS(Vec F, DM stokes_da, DM quadrature)
461: {
462: DM cda;
463: Vec coords;
464: const PetscScalar *_coords;
465: PetscInt u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
466: PetscInt p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
467: PetscInt nel, npe, eidx, i;
468: const PetscInt *element_list;
469: PetscScalar Fe[NODES_PER_EL * U_DOFS];
470: PetscScalar He[NODES_PER_EL * P_DOFS];
471: PetscScalar el_coords[NODES_PER_EL * NSD];
472: PetscScalar *q_rhs, *prop_fy;
473: Vec local_F;
474: PetscScalar *LA_F;
476: PetscFunctionBeginUser;
477: PetscCall(VecZeroEntries(F));
478: /* setup for coords */
479: PetscCall(DMGetCoordinateDM(stokes_da, &cda));
480: PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
481: PetscCall(VecGetArrayRead(coords, &_coords));
483: /* setup for coefficients */
484: PetscCall(DMSwarmGetField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
486: /* get access to the vector */
487: PetscCall(DMGetLocalVector(stokes_da, &local_F));
488: PetscCall(VecZeroEntries(local_F));
489: PetscCall(VecGetArray(local_F, &LA_F));
491: PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
492: for (eidx = 0; eidx < nel; eidx++) {
493: const PetscInt *element = &element_list[npe * eidx];
495: /* get coords for the element */
496: PetscCall(GetElementCoords(_coords, element, el_coords));
498: /* get coefficients for the element */
499: prop_fy = &q_rhs[GAUSS_POINTS * eidx];
501: /* initialise element stiffness matrix */
502: PetscCall(PetscMemzero(Fe, sizeof(Fe)));
503: PetscCall(PetscMemzero(He, sizeof(He)));
505: /* form element stiffness matrix */
506: LForm_MomentumRHS(Fe, el_coords, NULL, prop_fy);
508: /* insert element matrix into global matrix */
509: PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
511: for (i = 0; i < NODES_PER_EL * U_DOFS; i++) LA_F[u_eqn[i]] += Fe[i];
512: }
513: PetscCall(DMSwarmRestoreField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
514: PetscCall(VecRestoreArrayRead(coords, &_coords));
516: PetscCall(VecRestoreArray(local_F, &LA_F));
517: PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
518: PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
519: PetscCall(DMRestoreLocalVector(stokes_da, &local_F));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: PetscErrorCode DMSwarmPICInsertPointsCellwise(DM dm, DM dmc, PetscInt e, PetscInt npoints, PetscReal xi[], PetscBool proximity_initialization)
524: {
525: DMSwarmCellDM celldm;
526: PetscInt dim, nel, npe, q, k, d, ncurr, Nfc;
527: const PetscInt *element_list;
528: Vec coor;
529: const PetscScalar *_coor;
530: PetscReal **basis, *elcoor, *xp;
531: PetscReal *swarm_coor;
532: PetscInt *swarm_cellid;
533: const char **coordFields, *cellid;
535: PetscFunctionBeginUser;
536: PetscCall(DMGetDimension(dm, &dim));
537: PetscCall(DMDAGetElements(dmc, &nel, &npe, &element_list));
539: PetscCall(PetscMalloc1(dim * npoints, &xp));
540: PetscCall(PetscMalloc1(dim * npe, &elcoor));
541: PetscCall(PetscMalloc1(npoints, &basis));
542: for (q = 0; q < npoints; q++) {
543: PetscCall(PetscMalloc1(npe, &basis[q]));
545: switch (dim) {
546: case 1:
547: basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]);
548: basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]);
549: break;
550: case 2:
551: basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
552: basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
553: basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
554: basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
555: break;
557: case 3:
558: basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
559: basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
560: basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
561: basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
562: basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
563: basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
564: basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
565: basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
566: break;
567: }
568: }
570: PetscCall(DMGetCoordinatesLocal(dmc, &coor));
571: PetscCall(VecGetArrayRead(coor, &_coor));
572: /* compute and store the coordinates for the new points */
573: {
574: const PetscInt *element = &element_list[npe * e];
576: for (k = 0; k < npe; k++) {
577: for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]);
578: }
579: for (q = 0; q < npoints; q++) {
580: for (d = 0; d < dim; d++) xp[dim * q + d] = 0.0;
581: for (k = 0; k < npe; k++) {
582: for (d = 0; d < dim; d++) xp[dim * q + d] += basis[q][k] * elcoor[dim * k + d];
583: }
584: }
585: }
586: PetscCall(VecRestoreArrayRead(coor, &_coor));
587: PetscCall(DMDARestoreElements(dmc, &nel, &npe, &element_list));
589: PetscCall(DMSwarmGetLocalSize(dm, &ncurr));
590: PetscCall(DMSwarmAddNPoints(dm, npoints));
591: PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
592: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
593: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
594: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
596: if (proximity_initialization) {
597: PetscInt *nnlist;
598: PetscReal *coor_q, *coor_qn;
599: PetscInt npoints_e, *plist_e;
601: PetscCall(DMSwarmSortGetPointsPerCell(dm, e, &npoints_e, &plist_e));
603: PetscCall(PetscMalloc1(npoints, &nnlist));
604: /* find nearest neighbour points in this cell */
605: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
606: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
607: for (q = 0; q < npoints; q++) {
608: PetscInt qn, nearest_neighbour = -1;
609: PetscReal sep, min_sep = PETSC_MAX_REAL;
611: coor_q = &xp[dim * q];
612: for (qn = 0; qn < npoints_e; qn++) {
613: coor_qn = &swarm_coor[dim * plist_e[qn]];
614: sep = 0.0;
615: for (d = 0; d < dim; d++) sep += (coor_q[d] - coor_qn[d]) * (coor_q[d] - coor_qn[d]);
616: if (sep < min_sep) {
617: nearest_neighbour = plist_e[qn];
618: min_sep = sep;
619: }
620: }
621: PetscCheck(nearest_neighbour != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell %" PetscInt_FMT " is empty - cannot initialize using nearest neighbours", e);
622: nnlist[q] = nearest_neighbour;
623: }
624: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
625: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
627: /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
628: for (q = 0; q < npoints; q++) PetscCall(DMSwarmCopyPoint(dm, nnlist[q], ncurr + q));
629: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
630: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
631: for (q = 0; q < npoints; q++) {
632: /* set the coordinates */
633: for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
634: /* set the cell index */
635: swarm_cellid[ncurr + q] = e;
636: }
637: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
638: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
640: PetscCall(DMSwarmSortRestorePointsPerCell(dm, e, &npoints_e, &plist_e));
641: PetscCall(PetscFree(nnlist));
642: } else {
643: PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
644: PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
645: for (q = 0; q < npoints; q++) {
646: /* set the coordinates */
647: for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
648: /* set the cell index */
649: swarm_cellid[ncurr + q] = e;
650: }
651: PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&swarm_cellid));
652: PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&swarm_coor));
653: }
655: PetscCall(PetscFree(xp));
656: PetscCall(PetscFree(elcoor));
657: for (q = 0; q < npoints; q++) PetscCall(PetscFree(basis[q]));
658: PetscCall(PetscFree(basis));
659: PetscFunctionReturn(PETSC_SUCCESS);
660: }
662: PetscErrorCode MaterialPoint_PopulateCell(DM dm_vp, DM dm_mpoint)
663: {
664: PetscInt _npe, _nel, e, nel;
665: const PetscInt *element;
666: DM dmc;
667: PetscQuadrature quadrature;
668: const PetscReal *xi;
669: PetscInt npoints_q, cnt, cnt_g;
671: PetscFunctionBeginUser;
672: PetscCall(DMDAGetElements(dm_vp, &_nel, &_npe, &element));
673: nel = _nel;
674: PetscCall(DMDARestoreElements(dm_vp, &_nel, &_npe, &element));
676: PetscCall(PetscDTGaussTensorQuadrature(2, 1, 4, -1.0, 1.0, &quadrature));
677: PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xi, NULL));
678: PetscCall(DMSwarmGetCellDM(dm_mpoint, &dmc));
680: PetscCall(DMSwarmSortGetAccess(dm_mpoint));
682: cnt = 0;
683: for (e = 0; e < nel; e++) {
684: PetscInt npoints_per_cell;
686: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm_mpoint, e, &npoints_per_cell));
688: if (npoints_per_cell < 12) {
689: PetscCall(DMSwarmPICInsertPointsCellwise(dm_mpoint, dm_vp, e, npoints_q, (PetscReal *)xi, PETSC_TRUE));
690: cnt++;
691: }
692: }
693: PetscCallMPI(MPIU_Allreduce(&cnt, &cnt_g, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
694: if (cnt_g > 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... ....pop cont: adjusted %" PetscInt_FMT " cells\n", cnt_g));
696: PetscCall(DMSwarmSortRestoreAccess(dm_mpoint));
697: PetscCall(PetscQuadratureDestroy(&quadrature));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: PetscErrorCode MaterialPoint_AdvectRK1(DM dm_vp, Vec vp, PetscReal dt, DM dm_mpoint)
702: {
703: DMSwarmCellDM celldm;
704: Vec vp_l, coor_l;
705: const PetscScalar *LA_vp;
706: PetscInt i, p, e, npoints, nel, npe, Nfc;
707: PetscInt *mpfield_cell;
708: PetscReal *mpfield_coor;
709: const PetscInt *element_list;
710: const PetscInt *element;
711: PetscScalar xi_p[NSD], Ni[NODES_PER_EL];
712: const PetscScalar *LA_coor;
713: PetscScalar dx[NSD];
714: const char **coordFields, *cellid;
716: PetscFunctionBeginUser;
717: PetscCall(DMGetCoordinatesLocal(dm_vp, &coor_l));
718: PetscCall(VecGetArrayRead(coor_l, &LA_coor));
720: PetscCall(DMGetLocalVector(dm_vp, &vp_l));
721: PetscCall(DMGlobalToLocalBegin(dm_vp, vp, INSERT_VALUES, vp_l));
722: PetscCall(DMGlobalToLocalEnd(dm_vp, vp, INSERT_VALUES, vp_l));
723: PetscCall(VecGetArrayRead(vp_l, &LA_vp));
725: PetscCall(DMDAGetElements(dm_vp, &nel, &npe, &element_list));
726: PetscCall(DMSwarmGetLocalSize(dm_mpoint, &npoints));
727: PetscCall(DMSwarmGetCellDMActive(dm_mpoint, &celldm));
728: PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
729: PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
730: PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm_mpoint), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
731: PetscCall(DMSwarmGetField(dm_mpoint, coordFields[0], NULL, NULL, (void **)&mpfield_coor));
732: PetscCall(DMSwarmGetField(dm_mpoint, cellid, NULL, NULL, (void **)&mpfield_cell));
733: for (p = 0; p < npoints; p++) {
734: PetscReal *coor_p;
735: PetscScalar vel_n[NSD * NODES_PER_EL], vel_p[NSD];
736: const PetscScalar *x0;
737: const PetscScalar *x2;
739: e = mpfield_cell[p];
740: coor_p = &mpfield_coor[NSD * p];
741: element = &element_list[NODES_PER_EL * e];
743: /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
744: x0 = &LA_coor[NSD * element[0]];
745: x2 = &LA_coor[NSD * element[2]];
747: dx[0] = x2[0] - x0[0];
748: dx[1] = x2[1] - x0[1];
750: xi_p[0] = 2.0 * (coor_p[0] - x0[0]) / dx[0] - 1.0;
751: xi_p[1] = 2.0 * (coor_p[1] - x0[1]) / dx[1] - 1.0;
752: 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);
753: 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);
754: 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);
755: 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);
757: /* evaluate basis functions */
758: EvaluateBasis_Q1(xi_p, Ni);
760: /* get cell nodal velocities */
761: for (i = 0; i < NODES_PER_EL; i++) {
762: PetscInt nid;
764: nid = element[i];
765: vel_n[NSD * i + 0] = LA_vp[(NSD + 1) * nid + 0];
766: vel_n[NSD * i + 1] = LA_vp[(NSD + 1) * nid + 1];
767: }
769: /* interpolate velocity */
770: vel_p[0] = vel_p[1] = 0.0;
771: for (i = 0; i < NODES_PER_EL; i++) {
772: vel_p[0] += Ni[i] * vel_n[NSD * i + 0];
773: vel_p[1] += Ni[i] * vel_n[NSD * i + 1];
774: }
776: coor_p[0] += dt * PetscRealPart(vel_p[0]);
777: coor_p[1] += dt * PetscRealPart(vel_p[1]);
778: }
780: PetscCall(DMSwarmRestoreField(dm_mpoint, cellid, NULL, NULL, (void **)&mpfield_cell));
781: PetscCall(DMSwarmRestoreField(dm_mpoint, coordFields[0], NULL, NULL, (void **)&mpfield_coor));
782: PetscCall(DMDARestoreElements(dm_vp, &nel, &npe, &element_list));
783: PetscCall(VecRestoreArrayRead(vp_l, &LA_vp));
784: PetscCall(DMRestoreLocalVector(dm_vp, &vp_l));
785: PetscCall(VecRestoreArrayRead(coor_l, &LA_coor));
786: PetscFunctionReturn(PETSC_SUCCESS);
787: }
789: PetscErrorCode MaterialPoint_Interpolate(DM dm, Vec eta_v, Vec rho_v, DM dm_quadrature)
790: {
791: Vec eta_l, rho_l;
792: PetscScalar *_eta_l, *_rho_l;
793: PetscInt nqp, npe, nel;
794: PetscScalar qp_xi[GAUSS_POINTS][NSD];
795: PetscScalar qp_weight[GAUSS_POINTS];
796: PetscInt q, k, e;
797: PetscScalar Ni[GAUSS_POINTS][NODES_PER_EL];
798: const PetscInt *element_list;
799: PetscReal *q_eta, *q_rhs;
801: PetscFunctionBeginUser;
802: /* define quadrature rule */
803: CreateGaussQuadrature(&nqp, qp_xi, qp_weight);
804: for (q = 0; q < nqp; q++) EvaluateBasis_Q1(qp_xi[q], Ni[q]);
806: PetscCall(DMGetLocalVector(dm, &eta_l));
807: PetscCall(DMGetLocalVector(dm, &rho_l));
809: PetscCall(DMGlobalToLocalBegin(dm, eta_v, INSERT_VALUES, eta_l));
810: PetscCall(DMGlobalToLocalEnd(dm, eta_v, INSERT_VALUES, eta_l));
811: PetscCall(DMGlobalToLocalBegin(dm, rho_v, INSERT_VALUES, rho_l));
812: PetscCall(DMGlobalToLocalEnd(dm, rho_v, INSERT_VALUES, rho_l));
814: PetscCall(VecGetArray(eta_l, &_eta_l));
815: PetscCall(VecGetArray(rho_l, &_rho_l));
817: PetscCall(DMSwarmGetField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
818: PetscCall(DMSwarmGetField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
820: PetscCall(DMDAGetElements(dm, &nel, &npe, &element_list));
821: for (e = 0; e < nel; e++) {
822: PetscScalar eta_field_e[NODES_PER_EL];
823: PetscScalar rho_field_e[NODES_PER_EL];
824: const PetscInt *element = &element_list[4 * e];
826: for (k = 0; k < NODES_PER_EL; k++) {
827: eta_field_e[k] = _eta_l[element[k]];
828: rho_field_e[k] = _rho_l[element[k]];
829: }
831: for (q = 0; q < nqp; q++) {
832: PetscScalar eta_q, rho_q;
834: eta_q = rho_q = 0.0;
835: for (k = 0; k < NODES_PER_EL; k++) {
836: eta_q += Ni[q][k] * eta_field_e[k];
837: rho_q += Ni[q][k] * rho_field_e[k];
838: }
840: q_eta[nqp * e + q] = PetscRealPart(eta_q);
841: q_rhs[nqp * e + q] = PetscRealPart(rho_q);
842: }
843: }
844: PetscCall(DMDARestoreElements(dm, &nel, &npe, &element_list));
846: PetscCall(DMSwarmRestoreField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
847: PetscCall(DMSwarmRestoreField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
849: PetscCall(VecRestoreArray(rho_l, &_rho_l));
850: PetscCall(VecRestoreArray(eta_l, &_eta_l));
851: PetscCall(DMRestoreLocalVector(dm, &rho_l));
852: PetscCall(DMRestoreLocalVector(dm, &eta_l));
853: PetscFunctionReturn(PETSC_SUCCESS);
854: }
856: static PetscErrorCode SolveTimeDepStokes(PetscInt mx, PetscInt my)
857: {
858: DM dm_stokes, dm_coeff;
859: PetscInt u_dof, p_dof, dof, stencil_width;
860: Mat A, B;
861: PetscInt nel_local;
862: Vec eta_v, rho_v;
863: Vec f, X;
864: KSP ksp;
865: PC pc;
866: char filename[PETSC_MAX_PATH_LEN];
867: DM dms_quadrature, dms_mpoint;
868: PetscInt nel, npe, npoints;
869: const PetscInt *element_list;
870: PetscInt tk, nt, dump_freq;
871: PetscReal dt, dt_max = 0.0;
872: PetscReal vx[2], vy[2], max_v = 0.0, max_v_step, dh;
873: const char *fieldnames[] = {"eta", "rho"};
874: Vec pfields[2];
875: PetscInt ppcell = 1;
876: PetscReal time, delta_eta = 1.0;
877: PetscBool randomize_coords = PETSC_FALSE;
878: PetscReal randomize_fac = 0.25;
879: PetscBool no_view = PETSC_FALSE;
880: PetscBool isbddc;
882: PetscFunctionBeginUser;
883: /*
884: Generate the DMDA for the velocity and pressure spaces.
885: We use Q1 elements for both fields.
886: The Q1 FE basis on a regular mesh has a 9-point stencil (DMDA_STENCIL_BOX)
887: The number of nodes in each direction is mx+1, my+1
888: */
889: u_dof = U_DOFS; /* Vx, Vy - velocities */
890: p_dof = P_DOFS; /* p - pressure */
891: dof = u_dof + p_dof;
892: stencil_width = 1;
893: 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));
894: PetscCall(DMDASetElementType(dm_stokes, DMDA_ELEMENT_Q1));
895: PetscCall(DMSetMatType(dm_stokes, MATAIJ));
896: PetscCall(DMSetFromOptions(dm_stokes));
897: PetscCall(DMSetUp(dm_stokes));
898: PetscCall(DMDASetFieldName(dm_stokes, 0, "ux"));
899: PetscCall(DMDASetFieldName(dm_stokes, 1, "uy"));
900: PetscCall(DMDASetFieldName(dm_stokes, 2, "p"));
902: /* unit box [0,0.9142] x [0,1] */
903: PetscCall(DMDASetUniformCoordinates(dm_stokes, 0.0, 0.9142, 0.0, 1.0, 0., 0.));
904: dh = 1.0 / (PetscReal)mx;
906: /* Get local number of elements */
907: {
908: PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));
910: nel_local = nel;
912: PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
913: }
915: /* Create DMDA for representing scalar fields */
916: PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 1, &dm_coeff));
918: /* Create the swarm for storing quadrature point values */
919: PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_quadrature));
920: PetscCall(DMSetType(dms_quadrature, DMSWARM));
921: PetscCall(DMSetDimension(dms_quadrature, 2));
922: PetscCall(PetscObjectSetName((PetscObject)dms_quadrature, "Quadrature Swarm"));
924: /* Register fields for viscosity and density on the quadrature points */
925: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "eta_q", 1, PETSC_REAL));
926: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "rho_q", 1, PETSC_REAL));
927: PetscCall(DMSwarmFinalizeFieldRegister(dms_quadrature));
928: PetscCall(DMSwarmSetLocalSizes(dms_quadrature, nel_local * GAUSS_POINTS, 0));
930: /* Create the material point swarm */
931: PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_mpoint));
932: PetscCall(DMSetType(dms_mpoint, DMSWARM));
933: PetscCall(DMSetDimension(dms_mpoint, 2));
934: PetscCall(PetscObjectSetName((PetscObject)dms_mpoint, "Material Point Swarm"));
936: /* Configure the material point swarm to be of type Particle-In-Cell */
937: PetscCall(DMSwarmSetType(dms_mpoint, DMSWARM_PIC));
939: /*
940: Specify the DM to use for point location and projections
941: within the context of a PIC scheme
942: */
943: PetscCall(DMSwarmSetCellDM(dms_mpoint, dm_coeff));
945: /* Register fields for viscosity and density */
946: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "eta", 1, PETSC_REAL));
947: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "rho", 1, PETSC_REAL));
948: PetscCall(DMSwarmFinalizeFieldRegister(dms_mpoint));
950: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL));
951: PetscCall(DMSwarmSetLocalSizes(dms_mpoint, nel_local * ppcell, 100));
953: /*
954: Layout the material points in space using the cell DM.
955: Particle coordinates are defined by cell wise using different methods.
956: - DMSWARMPIC_LAYOUT_GAUSS defines particles coordinates at the positions
957: corresponding to a Gauss quadrature rule with
958: ppcell points in each direction.
959: - DMSWARMPIC_LAYOUT_REGULAR defines particle coordinates at the centoid of
960: ppcell x ppcell quadralaterals defined within the
961: reference element.
962: - DMSWARMPIC_LAYOUT_SUBDIVISION defines particles coordinates at the centroid
963: of each quadralateral obtained by sub-dividing
964: the reference element cell ppcell times.
965: */
966: PetscCall(DMSwarmInsertPointsUsingCellDM(dms_mpoint, DMSWARMPIC_LAYOUT_SUBDIVISION, ppcell));
968: /*
969: Defne a high resolution layer of material points across the material interface
970: */
971: {
972: PetscInt npoints_dir_x[2];
973: PetscReal min[2], max[2];
975: npoints_dir_x[0] = (PetscInt)(0.9142 / (0.05 * dh));
976: npoints_dir_x[1] = (PetscInt)((0.25 - 0.15) / (0.05 * dh));
977: min[0] = 0.0;
978: max[0] = 0.9142;
979: min[1] = 0.05;
980: max[1] = 0.35;
981: PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
982: }
984: /*
985: Define a high resolution layer of material points near the surface of the domain
986: to deal with weakly compressible Q1-Q1 elements. These elements "self compact"
987: when applied to buoyancy driven flow. The error in div(u) is O(h).
988: */
989: {
990: PetscInt npoints_dir_x[2];
991: PetscReal min[2], max[2];
993: npoints_dir_x[0] = (PetscInt)(0.9142 / (0.25 * dh));
994: npoints_dir_x[1] = (PetscInt)(3.0 * dh / (0.25 * dh));
995: min[0] = 0.0;
996: max[0] = 0.9142;
997: min[1] = 1.0 - 3.0 * dh;
998: max[1] = 1.0 - 0.0001;
999: PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
1000: }
1002: PetscCall(DMView(dms_mpoint, PETSC_VIEWER_STDOUT_WORLD));
1004: /* Define initial material properties on each particle in the material point swarm */
1005: PetscCall(PetscOptionsGetReal(NULL, NULL, "-delta_eta", &delta_eta, NULL));
1006: PetscCall(PetscOptionsGetBool(NULL, NULL, "-randomize_coords", &randomize_coords, NULL));
1007: PetscCall(PetscOptionsGetReal(NULL, NULL, "-randomize_fac", &randomize_fac, NULL));
1008: PetscCheck(randomize_fac <= 1.0, PETSC_COMM_WORLD, PETSC_ERR_USER, "The value of -randomize_fac should be <= 1.0");
1009: {
1010: PetscReal *array_x, *array_e, *array_r;
1011: PetscInt p;
1012: PetscRandom r;
1013: PetscMPIInt rank;
1015: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1017: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
1018: PetscCall(PetscRandomSetInterval(r, -randomize_fac * dh, randomize_fac * dh));
1019: PetscCall(PetscRandomSetSeed(r, (unsigned long)rank));
1020: PetscCall(PetscRandomSeed(r));
1022: PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));
1024: /*
1025: Fetch the registered data from the material point DMSwarm.
1026: The fields "eta" and "rho" were registered by this example.
1027: The field identified by the variable DMSwarmPICField_coor
1028: was registered by the DMSwarm implementation when the function
1029: DMSwarmSetType(dms_mpoint,DMSWARM_PIC)
1030: was called. The returned array defines the coordinates of each
1031: material point in the point swarm.
1032: */
1033: PetscCall(DMSwarmGetField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));
1034: PetscCall(DMSwarmGetField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
1035: PetscCall(DMSwarmGetField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));
1037: PetscCall(DMSwarmGetLocalSize(dms_mpoint, &npoints));
1038: for (p = 0; p < npoints; p++) {
1039: PetscReal x_p[2], rr[2];
1041: if (randomize_coords) {
1042: PetscCall(PetscRandomGetValueReal(r, &rr[0]));
1043: PetscCall(PetscRandomGetValueReal(r, &rr[1]));
1044: array_x[2 * p + 0] += rr[0];
1045: array_x[2 * p + 1] += rr[1];
1046: }
1048: /* Get the coordinates of point, p */
1049: x_p[0] = array_x[2 * p + 0];
1050: x_p[1] = array_x[2 * p + 1];
1052: if (x_p[1] < (0.2 + 0.02 * PetscCosReal(PETSC_PI * x_p[0] / 0.9142))) {
1053: /* Material properties below the interface */
1054: array_e[p] = 1.0 * (1.0 / delta_eta);
1055: array_r[p] = 0.0;
1056: } else {
1057: /* Material properties above the interface */
1058: array_e[p] = 1.0;
1059: array_r[p] = 1.0;
1060: }
1061: }
1063: /*
1064: Restore the fetched data fields from the material point DMSwarm.
1065: Calling the Restore function invalidates the points array_r, array_e, array_x
1066: by setting them to NULL.
1067: */
1068: PetscCall(DMSwarmRestoreField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));
1069: PetscCall(DMSwarmRestoreField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
1070: PetscCall(DMSwarmRestoreField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));
1072: PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
1073: PetscCall(PetscRandomDestroy(&r));
1074: }
1076: /*
1077: If the particle coordinates where randomly shifted, they may have crossed into another
1078: element, or into another sub-domain. To account for this we call the Migrate function.
1079: */
1080: if (randomize_coords) PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));
1082: PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL));
1083: if (!no_view) PetscCall(DMSwarmViewXDMF(dms_mpoint, "ic_coeff_dms.xmf"));
1085: /* project the swarm properties */
1086: PetscCall(DMCreateGlobalVector(dm_coeff, &pfields[0]));
1087: PetscCall(DMCreateGlobalVector(dm_coeff, &pfields[1]));
1088: PetscCall(DMSwarmProjectFields(dms_mpoint, NULL, 2, fieldnames, pfields, SCATTER_FORWARD));
1089: eta_v = pfields[0];
1090: rho_v = pfields[1];
1091: PetscCall(PetscObjectSetName((PetscObject)eta_v, "eta"));
1092: PetscCall(PetscObjectSetName((PetscObject)rho_v, "rho"));
1093: PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));
1095: /* view projected coefficients eta and rho */
1096: if (!no_view) {
1097: PetscViewer viewer;
1099: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1100: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1101: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1102: PetscCall(PetscViewerFileSetName(viewer, "ic_coeff_dmda.vts"));
1103: PetscCall(VecView(eta_v, viewer));
1104: PetscCall(VecView(rho_v, viewer));
1105: PetscCall(PetscViewerDestroy(&viewer));
1106: }
1108: PetscCall(DMCreateMatrix(dm_stokes, &A));
1109: PetscCall(DMCreateMatrix(dm_stokes, &B));
1110: PetscCall(DMCreateGlobalVector(dm_stokes, &f));
1111: PetscCall(DMCreateGlobalVector(dm_stokes, &X));
1113: PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
1114: PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
1115: PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));
1117: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
1118: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));
1120: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
1121: PetscCall(KSPSetOptionsPrefix(ksp, "stokes_"));
1122: PetscCall(KSPSetDM(ksp, dm_stokes));
1123: PetscCall(KSPSetDMActive(ksp, PETSC_FALSE));
1124: PetscCall(KSPSetOperators(ksp, A, B));
1125: PetscCall(KSPSetFromOptions(ksp));
1126: PetscCall(KSPGetPC(ksp, &pc));
1127: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &isbddc));
1128: if (isbddc) PetscCall(KSPSetOperators(ksp, A, A));
1130: /* Define u-v-p indices for fieldsplit */
1131: {
1132: PC pc;
1133: const PetscInt ufields[] = {0, 1}, pfields[1] = {2};
1135: PetscCall(KSPGetPC(ksp, &pc));
1136: PetscCall(PCFieldSplitSetBlockSize(pc, 3));
1137: PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
1138: PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
1139: }
1141: /* If using a fieldsplit preconditioner, attach a DMDA to the velocity split so that geometric multigrid can be used */
1142: {
1143: PC pc, pc_u;
1144: KSP *sub_ksp, ksp_u;
1145: PetscInt nsplits;
1146: DM dm_u;
1147: PetscBool is_pcfs;
1149: PetscCall(KSPGetPC(ksp, &pc));
1151: is_pcfs = PETSC_FALSE;
1152: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs));
1154: if (is_pcfs) {
1155: PetscCall(KSPSetUp(ksp));
1156: PetscCall(KSPGetPC(ksp, &pc));
1157: PetscCall(PCFieldSplitGetSubKSP(pc, &nsplits, &sub_ksp));
1158: ksp_u = sub_ksp[0];
1159: PetscCall(PetscFree(sub_ksp));
1161: if (nsplits == 2) {
1162: PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 2, &dm_u));
1164: PetscCall(KSPSetDM(ksp_u, dm_u));
1165: PetscCall(KSPSetDMActive(ksp_u, PETSC_FALSE));
1166: PetscCall(DMDestroy(&dm_u));
1168: /* enforce galerkin coarse grids be used */
1169: PetscCall(KSPGetPC(ksp_u, &pc_u));
1170: PetscCall(PCMGSetGalerkin(pc_u, PC_MG_GALERKIN_PMAT));
1171: }
1172: }
1173: }
1175: dump_freq = 10;
1176: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dump_freq", &dump_freq, NULL));
1177: nt = 10;
1178: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nt", &nt, NULL));
1179: time = 0.0;
1180: for (tk = 1; tk <= nt; tk++) {
1181: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... assemble\n"));
1182: PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
1183: PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
1184: PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));
1186: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... bc imposition\n"));
1187: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
1188: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));
1190: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... solve\n"));
1191: PetscCall(KSPSetOperators(ksp, A, isbddc ? A : B));
1192: PetscCall(KSPSolve(ksp, f, X));
1194: PetscCall(VecStrideMax(X, 0, NULL, &vx[1]));
1195: PetscCall(VecStrideMax(X, 1, NULL, &vy[1]));
1196: PetscCall(VecStrideMin(X, 0, NULL, &vx[0]));
1197: PetscCall(VecStrideMin(X, 1, NULL, &vy[0]));
1199: max_v_step = PetscMax(vx[0], vx[1]);
1200: max_v_step = PetscMax(max_v_step, vy[0]);
1201: max_v_step = PetscMax(max_v_step, vy[1]);
1202: max_v = PetscMax(max_v, max_v_step);
1204: dt_max = 2.0;
1205: dt = 0.5 * (dh / max_v_step);
1206: 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));
1207: dt = PetscMin(dt_max, dt);
1209: /* advect */
1210: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... advect\n"));
1211: PetscCall(MaterialPoint_AdvectRK1(dm_stokes, X, dt, dms_mpoint));
1213: /* migrate */
1214: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... migrate\n"));
1215: PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));
1217: /* update cell population */
1218: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... populate cells\n"));
1219: PetscCall(MaterialPoint_PopulateCell(dm_stokes, dms_mpoint));
1221: /* update coefficients on quadrature points */
1222: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... project\n"));
1223: PetscCall(DMSwarmProjectFields(dms_mpoint, NULL, 2, fieldnames, pfields, SCATTER_FORWARD));
1224: eta_v = pfields[0];
1225: rho_v = pfields[1];
1226: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... interp\n"));
1227: PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));
1229: if (tk % dump_freq == 0) {
1230: PetscViewer viewer;
1232: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... write XDMF, VTS\n"));
1233: PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_coeff_dms.xmf", tk));
1234: PetscCall(DMSwarmViewXDMF(dms_mpoint, filename));
1236: PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_vp_dm.vts", tk));
1237: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1238: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1239: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1240: PetscCall(PetscViewerFileSetName(viewer, filename));
1241: PetscCall(VecView(X, viewer));
1242: PetscCall(PetscViewerDestroy(&viewer));
1243: }
1244: time += dt;
1245: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step %" PetscInt_FMT " : time %1.2e \n", tk, (double)time));
1246: }
1248: PetscCall(KSPDestroy(&ksp));
1249: PetscCall(VecDestroy(&X));
1250: PetscCall(VecDestroy(&f));
1251: PetscCall(MatDestroy(&A));
1252: PetscCall(MatDestroy(&B));
1253: PetscCall(VecDestroy(&eta_v));
1254: PetscCall(VecDestroy(&rho_v));
1256: PetscCall(DMDestroy(&dms_mpoint));
1257: PetscCall(DMDestroy(&dms_quadrature));
1258: PetscCall(DMDestroy(&dm_coeff));
1259: PetscCall(DMDestroy(&dm_stokes));
1260: PetscFunctionReturn(PETSC_SUCCESS);
1261: }
1263: /*
1264: <sequential run>
1265: ./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
1266: */
1267: int main(int argc, char **args)
1268: {
1269: PetscInt mx, my;
1270: PetscBool set = PETSC_FALSE;
1272: PetscFunctionBeginUser;
1273: PetscCall(PetscInitialize(&argc, &args, NULL, help));
1274: mx = my = 10;
1275: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1276: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1277: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mxy", &mx, &set));
1278: if (set) my = mx;
1279: PetscCall(SolveTimeDepStokes(mx, my));
1280: PetscCall(PetscFinalize());
1281: return 0;
1282: }
1284: /* -------------------------- helpers for boundary conditions -------------------------------- */
1285: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
1286: {
1287: DM cda;
1288: Vec coords;
1289: PetscInt si, sj, nx, ny, i, j;
1290: PetscInt M, N;
1291: DMDACoor2d **_coords;
1292: const PetscInt *g_idx;
1293: PetscInt *bc_global_ids;
1294: PetscScalar *bc_vals;
1295: PetscInt nbcs;
1296: PetscInt n_dofs;
1297: ISLocalToGlobalMapping ltogm;
1299: PetscFunctionBeginUser;
1300: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1301: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1303: PetscCall(DMGetCoordinateDM(da, &cda));
1304: PetscCall(DMGetCoordinatesLocal(da, &coords));
1305: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1306: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1307: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1309: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1310: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1312: /* init the entries to -1 so VecSetValues will ignore them */
1313: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1315: i = nx - 1;
1316: for (j = 0; j < ny; j++) {
1317: PetscInt local_id;
1319: local_id = i + j * nx;
1321: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1323: bc_vals[j] = 0.0;
1324: }
1325: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1326: nbcs = 0;
1327: if ((si + nx) == (M)) nbcs = ny;
1329: if (b) {
1330: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1331: PetscCall(VecAssemblyBegin(b));
1332: PetscCall(VecAssemblyEnd(b));
1333: }
1334: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1336: PetscCall(PetscFree(bc_vals));
1337: PetscCall(PetscFree(bc_global_ids));
1339: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1340: PetscFunctionReturn(PETSC_SUCCESS);
1341: }
1343: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
1344: {
1345: DM cda;
1346: Vec coords;
1347: PetscInt si, sj, nx, ny, i, j;
1348: PetscInt M, N;
1349: DMDACoor2d **_coords;
1350: const PetscInt *g_idx;
1351: PetscInt *bc_global_ids;
1352: PetscScalar *bc_vals;
1353: PetscInt nbcs;
1354: PetscInt n_dofs;
1355: ISLocalToGlobalMapping ltogm;
1357: PetscFunctionBeginUser;
1358: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1359: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1361: PetscCall(DMGetCoordinateDM(da, &cda));
1362: PetscCall(DMGetCoordinatesLocal(da, &coords));
1363: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1364: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1365: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1367: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1368: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1370: /* init the entries to -1 so VecSetValues will ignore them */
1371: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1373: i = 0;
1374: for (j = 0; j < ny; j++) {
1375: PetscInt local_id;
1377: local_id = i + j * nx;
1379: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1381: bc_vals[j] = 0.0;
1382: }
1383: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1384: nbcs = 0;
1385: if (si == 0) nbcs = ny;
1387: if (b) {
1388: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1389: PetscCall(VecAssemblyBegin(b));
1390: PetscCall(VecAssemblyEnd(b));
1391: }
1393: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1395: PetscCall(PetscFree(bc_vals));
1396: PetscCall(PetscFree(bc_global_ids));
1398: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1399: PetscFunctionReturn(PETSC_SUCCESS);
1400: }
1402: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
1403: {
1404: DM cda;
1405: Vec coords;
1406: PetscInt si, sj, nx, ny, i, j;
1407: PetscInt M, N;
1408: DMDACoor2d **_coords;
1409: const PetscInt *g_idx;
1410: PetscInt *bc_global_ids;
1411: PetscScalar *bc_vals;
1412: PetscInt nbcs;
1413: PetscInt n_dofs;
1414: ISLocalToGlobalMapping ltogm;
1416: PetscFunctionBeginUser;
1417: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1418: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1420: PetscCall(DMGetCoordinateDM(da, &cda));
1421: PetscCall(DMGetCoordinatesLocal(da, &coords));
1422: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1423: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1424: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1426: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1427: PetscCall(PetscMalloc1(nx, &bc_vals));
1429: /* init the entries to -1 so VecSetValues will ignore them */
1430: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1432: j = ny - 1;
1433: for (i = 0; i < nx; i++) {
1434: PetscInt local_id;
1436: local_id = i + j * nx;
1438: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1440: bc_vals[i] = 0.0;
1441: }
1442: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1443: nbcs = 0;
1444: if ((sj + ny) == (N)) nbcs = nx;
1446: if (b) {
1447: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1448: PetscCall(VecAssemblyBegin(b));
1449: PetscCall(VecAssemblyEnd(b));
1450: }
1451: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, NULL, NULL));
1453: PetscCall(PetscFree(bc_vals));
1454: PetscCall(PetscFree(bc_global_ids));
1456: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1457: PetscFunctionReturn(PETSC_SUCCESS);
1458: }
1460: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
1461: {
1462: DM cda;
1463: Vec coords;
1464: PetscInt si, sj, nx, ny, i, j;
1465: PetscInt M, N;
1466: DMDACoor2d **_coords;
1467: const PetscInt *g_idx;
1468: PetscInt *bc_global_ids;
1469: PetscScalar *bc_vals;
1470: PetscInt nbcs;
1471: PetscInt n_dofs;
1472: ISLocalToGlobalMapping ltogm;
1474: PetscFunctionBeginUser;
1475: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1476: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1478: PetscCall(DMGetCoordinateDM(da, &cda));
1479: PetscCall(DMGetCoordinatesLocal(da, &coords));
1480: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1481: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1482: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1484: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1485: PetscCall(PetscMalloc1(nx, &bc_vals));
1487: /* init the entries to -1 so VecSetValues will ignore them */
1488: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1490: j = 0;
1491: for (i = 0; i < nx; i++) {
1492: PetscInt local_id;
1494: local_id = i + j * nx;
1496: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1498: bc_vals[i] = 0.0;
1499: }
1500: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1501: nbcs = 0;
1502: if (sj == 0) nbcs = nx;
1504: if (b) {
1505: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1506: PetscCall(VecAssemblyBegin(b));
1507: PetscCall(VecAssemblyEnd(b));
1508: }
1509: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1511: PetscCall(PetscFree(bc_vals));
1512: PetscCall(PetscFree(bc_global_ids));
1514: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1515: PetscFunctionReturn(PETSC_SUCCESS);
1516: }
1518: /*
1519: Impose free slip boundary conditions on the left/right faces: u_i n_i = 0, tau_{ij} t_j = 0
1520: Impose no slip boundray conditions on the top/bottom faces: u_i n_i = 0, u_i t_i = 0
1521: */
1522: static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes, Mat A, Vec f)
1523: {
1524: PetscFunctionBeginUser;
1525: PetscCall(BCApplyZero_NORTH(dm_stokes, 0, A, f));
1526: PetscCall(BCApplyZero_NORTH(dm_stokes, 1, A, f));
1527: PetscCall(BCApplyZero_EAST(dm_stokes, 0, A, f));
1528: PetscCall(BCApplyZero_SOUTH(dm_stokes, 0, A, f));
1529: PetscCall(BCApplyZero_SOUTH(dm_stokes, 1, A, f));
1530: PetscCall(BCApplyZero_WEST(dm_stokes, 0, A, f));
1531: PetscFunctionReturn(PETSC_SUCCESS);
1532: }
1534: /*TEST
1536: test:
1537: suffix: 1
1538: args: -no_view
1539: requires: !complex double
1540: filter: grep -v atomic
1541: filter_output: grep -v atomic
1542: test:
1543: suffix: 1_matis
1544: requires: !complex double
1545: args: -no_view -dm_mat_type is
1546: filter: grep -v atomic
1547: filter_output: grep -v atomic
1548: testset:
1549: nsize: 4
1550: requires: !complex double
1551: 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
1552: filter: grep -v atomic
1553: filter_output: grep -v atomic
1554: test:
1555: suffix: fetidp
1556: args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
1557: test:
1558: suffix: fetidp_lumped
1559: 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
1560: test:
1561: suffix: fetidp_saddlepoint
1562: 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
1563: test:
1564: suffix: fetidp_saddlepoint_lumped
1565: 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
1566: TEST*/