Actual source code: ex70.c
1: static char help[] = "------------------------------------------------------------------------------------------------------------------------------ \n\
2: Solves the time-dependent incompressible, variable viscosity Stokes equation in 2D driven by buouyancy variations. \n\
3: Time-dependence is introduced by evolving the density (rho) and viscosity (eta) according to \n\
4: D \\rho / Dt = 0 and D \\eta / Dt = 0 \n\
5: The Stokes problem is discretized using Q1-Q1 finite elements, stabilized with Bochev's polynomial projection method. \n\
6: The hyperbolic evolution equation for density is discretized using a variant of the Particle-In-Cell (PIC) method. \n\
7: The DMDA object is used to define the FE problem, whilst DMSwarm provides support for the PIC method. \n\
8: Material points (particles) store density and viscosity. The particles are advected with the fluid velocity using RK1. \n\
9: At each time step, the value of density and viscosity stored on each particle are projected into a Q1 function space \n\
10: and then interpolated onto the Gauss quadrature points. \n\
11: The model problem defined in this example is the iso-viscous Rayleigh-Taylor instability (case 1a) from: \n\
12: \"A comparison of methods for the modeling of thermochemical convection\" \n\
13: P.E. van Keken, S.D. King, H. Schmeling, U.R. Christensen, D. Neumeister and M.-P. Doin, \n\
14: Journal of Geophysical Research, vol 102 (B10), 477--499 (1997) \n\
15: Note that whilst the model problem defined is for an iso-viscous, the implementation in this example supports \n\
16: variable viscoity formulations. \n\
17: This example is based on src/ksp/ksp/tutorials/ex43.c \n\
18: Options: \n\
19: -mx : Number of elements in the x-direction \n\
20: -my : Number of elements in the y-direction \n\
21: -mxy : Number of elements in the x- and y-directions \n\
22: -nt : Number of time steps \n\
23: -dump_freq : Frequency of output file creation \n\
24: -ppcell : Number of times the reference cell is sub-divided \n\
25: -randomize_coords : Apply a random shift to each particle coordinate in the range [-fac*dh,0.fac*dh] \n\
26: -randomize_fac : Set the scaling factor for the random shift (default = 0.25)\n";
28: /* Contributed by Dave May */
30: #include <petscksp.h>
31: #include <petscdm.h>
32: #include <petscdmda.h>
33: #include <petscdmswarm.h>
34: #include <petsc/private/dmimpl.h>
36: static PetscErrorCode DMDAApplyBoundaryConditions(DM, Mat, Vec);
38: #define NSD 2 /* number of spatial dimensions */
39: #define NODES_PER_EL 4 /* nodes per element */
40: #define U_DOFS 2 /* degrees of freedom per velocity node */
41: #define P_DOFS 1 /* degrees of freedom per pressure node */
42: #define GAUSS_POINTS 4
44: static void EvaluateBasis_Q1(PetscScalar _xi[], PetscScalar N[])
45: {
46: PetscScalar xi = _xi[0];
47: PetscScalar eta = _xi[1];
49: N[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
50: N[1] = 0.25 * (1.0 + xi) * (1.0 - eta);
51: N[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
52: N[3] = 0.25 * (1.0 - xi) * (1.0 + eta);
53: }
55: static void EvaluateBasisDerivatives_Q1(PetscScalar _xi[], PetscScalar dN[][NODES_PER_EL])
56: {
57: PetscScalar xi = _xi[0];
58: PetscScalar eta = _xi[1];
60: dN[0][0] = -0.25 * (1.0 - eta);
61: dN[0][1] = 0.25 * (1.0 - eta);
62: dN[0][2] = 0.25 * (1.0 + eta);
63: dN[0][3] = -0.25 * (1.0 + eta);
65: dN[1][0] = -0.25 * (1.0 - xi);
66: dN[1][1] = -0.25 * (1.0 + xi);
67: dN[1][2] = 0.25 * (1.0 + xi);
68: dN[1][3] = 0.25 * (1.0 - xi);
69: }
71: static void EvaluateDerivatives(PetscScalar dN[][NODES_PER_EL], PetscScalar dNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
72: {
73: PetscScalar J00, J01, J10, J11, J;
74: PetscScalar iJ00, iJ01, iJ10, iJ11;
75: PetscInt i;
77: J00 = J01 = J10 = J11 = 0.0;
78: for (i = 0; i < NODES_PER_EL; i++) {
79: PetscScalar cx = coords[2 * i];
80: PetscScalar cy = coords[2 * i + 1];
82: J00 += dN[0][i] * cx; /* J_xx = dx/dxi */
83: J01 += dN[0][i] * cy; /* J_xy = dy/dxi */
84: J10 += dN[1][i] * cx; /* J_yx = dx/deta */
85: J11 += dN[1][i] * cy; /* J_yy = dy/deta */
86: }
87: J = (J00 * J11) - (J01 * J10);
89: iJ00 = J11 / J;
90: iJ01 = -J01 / J;
91: iJ10 = -J10 / J;
92: iJ11 = J00 / J;
94: for (i = 0; i < NODES_PER_EL; i++) {
95: dNx[0][i] = dN[0][i] * iJ00 + dN[1][i] * iJ01;
96: dNx[1][i] = dN[0][i] * iJ10 + dN[1][i] * iJ11;
97: }
99: if (det_J) *det_J = J;
100: }
102: static void CreateGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
103: {
104: *ngp = 4;
105: gp_xi[0][0] = -0.57735026919;
106: gp_xi[0][1] = -0.57735026919;
107: gp_xi[1][0] = -0.57735026919;
108: gp_xi[1][1] = 0.57735026919;
109: gp_xi[2][0] = 0.57735026919;
110: gp_xi[2][1] = 0.57735026919;
111: gp_xi[3][0] = 0.57735026919;
112: gp_xi[3][1] = -0.57735026919;
113: gp_weight[0] = 1.0;
114: gp_weight[1] = 1.0;
115: gp_weight[2] = 1.0;
116: gp_weight[3] = 1.0;
117: }
119: static PetscErrorCode DMDAGetElementEqnums_up(const PetscInt element[], PetscInt s_u[], PetscInt s_p[])
120: {
121: 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: PetscInt dim, nel, npe, q, k, d, ncurr;
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;
533: PetscFunctionBeginUser;
534: PetscCall(DMGetDimension(dm, &dim));
535: PetscCall(DMDAGetElements(dmc, &nel, &npe, &element_list));
537: PetscCall(PetscMalloc1(dim * npoints, &xp));
538: PetscCall(PetscMalloc1(dim * npe, &elcoor));
539: PetscCall(PetscMalloc1(npoints, &basis));
540: for (q = 0; q < npoints; q++) {
541: PetscCall(PetscMalloc1(npe, &basis[q]));
543: switch (dim) {
544: case 1:
545: basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]);
546: basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]);
547: break;
548: case 2:
549: basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
550: basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
551: basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
552: basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
553: break;
555: case 3:
556: basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
557: basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
558: basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
559: basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
560: basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
561: basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
562: basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
563: basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
564: break;
565: }
566: }
568: PetscCall(DMGetCoordinatesLocal(dmc, &coor));
569: PetscCall(VecGetArrayRead(coor, &_coor));
570: /* compute and store the coordinates for the new points */
571: {
572: const PetscInt *element = &element_list[npe * e];
574: for (k = 0; k < npe; k++) {
575: for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]);
576: }
577: for (q = 0; q < npoints; q++) {
578: for (d = 0; d < dim; d++) xp[dim * q + d] = 0.0;
579: for (k = 0; k < npe; k++) {
580: for (d = 0; d < dim; d++) xp[dim * q + d] += basis[q][k] * elcoor[dim * k + d];
581: }
582: }
583: }
584: PetscCall(VecRestoreArrayRead(coor, &_coor));
585: PetscCall(DMDARestoreElements(dmc, &nel, &npe, &element_list));
587: PetscCall(DMSwarmGetLocalSize(dm, &ncurr));
588: PetscCall(DMSwarmAddNPoints(dm, npoints));
590: if (proximity_initialization) {
591: PetscInt *nnlist;
592: PetscReal *coor_q, *coor_qn;
593: PetscInt npoints_e, *plist_e;
595: PetscCall(DMSwarmSortGetPointsPerCell(dm, e, &npoints_e, &plist_e));
597: PetscCall(PetscMalloc1(npoints, &nnlist));
598: /* find nearest neighbour points in this cell */
599: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
600: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
601: for (q = 0; q < npoints; q++) {
602: PetscInt qn, nearest_neighbour = -1;
603: PetscReal sep, min_sep = PETSC_MAX_REAL;
605: coor_q = &xp[dim * q];
606: for (qn = 0; qn < npoints_e; qn++) {
607: coor_qn = &swarm_coor[dim * plist_e[qn]];
608: sep = 0.0;
609: for (d = 0; d < dim; d++) sep += (coor_q[d] - coor_qn[d]) * (coor_q[d] - coor_qn[d]);
610: if (sep < min_sep) {
611: nearest_neighbour = plist_e[qn];
612: min_sep = sep;
613: }
614: }
615: PetscCheck(nearest_neighbour != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell %" PetscInt_FMT " is empty - cannot initialize using nearest neighbours", e);
616: nnlist[q] = nearest_neighbour;
617: }
618: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
619: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
621: /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
622: for (q = 0; q < npoints; q++) PetscCall(DMSwarmCopyPoint(dm, nnlist[q], ncurr + q));
623: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
624: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
625: for (q = 0; q < npoints; q++) {
626: /* set the coordinates */
627: for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
628: /* set the cell index */
629: swarm_cellid[ncurr + q] = e;
630: }
631: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
632: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
634: PetscCall(PetscFree(plist_e));
635: PetscCall(PetscFree(nnlist));
636: } else {
637: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
638: PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
639: for (q = 0; q < npoints; q++) {
640: /* set the coordinates */
641: for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
642: /* set the cell index */
643: swarm_cellid[ncurr + q] = e;
644: }
645: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
646: PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
647: }
649: PetscCall(PetscFree(xp));
650: PetscCall(PetscFree(elcoor));
651: for (q = 0; q < npoints; q++) PetscCall(PetscFree(basis[q]));
652: PetscCall(PetscFree(basis));
653: PetscFunctionReturn(PETSC_SUCCESS);
654: }
656: PetscErrorCode MaterialPoint_PopulateCell(DM dm_vp, DM dm_mpoint)
657: {
658: PetscInt _npe, _nel, e, nel;
659: const PetscInt *element;
660: DM dmc;
661: PetscQuadrature quadrature;
662: const PetscReal *xi;
663: PetscInt npoints_q, cnt, cnt_g;
665: PetscFunctionBeginUser;
666: PetscCall(DMDAGetElements(dm_vp, &_nel, &_npe, &element));
667: nel = _nel;
668: PetscCall(DMDARestoreElements(dm_vp, &_nel, &_npe, &element));
670: PetscCall(PetscDTGaussTensorQuadrature(2, 1, 4, -1.0, 1.0, &quadrature));
671: PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xi, NULL));
672: PetscCall(DMSwarmGetCellDM(dm_mpoint, &dmc));
674: PetscCall(DMSwarmSortGetAccess(dm_mpoint));
676: cnt = 0;
677: for (e = 0; e < nel; e++) {
678: PetscInt npoints_per_cell;
680: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm_mpoint, e, &npoints_per_cell));
682: if (npoints_per_cell < 12) {
683: PetscCall(DMSwarmPICInsertPointsCellwise(dm_mpoint, dm_vp, e, npoints_q, (PetscReal *)xi, PETSC_TRUE));
684: cnt++;
685: }
686: }
687: PetscCall(MPIU_Allreduce(&cnt, &cnt_g, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
688: if (cnt_g > 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... ....pop cont: adjusted %" PetscInt_FMT " cells\n", cnt_g));
690: PetscCall(DMSwarmSortRestoreAccess(dm_mpoint));
691: PetscCall(PetscQuadratureDestroy(&quadrature));
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: PetscErrorCode MaterialPoint_AdvectRK1(DM dm_vp, Vec vp, PetscReal dt, DM dm_mpoint)
696: {
697: Vec vp_l, coor_l;
698: const PetscScalar *LA_vp;
699: PetscInt i, p, e, npoints, nel, npe;
700: PetscInt *mpfield_cell;
701: PetscReal *mpfield_coor;
702: const PetscInt *element_list;
703: const PetscInt *element;
704: PetscScalar xi_p[NSD], Ni[NODES_PER_EL];
705: const PetscScalar *LA_coor;
706: PetscScalar dx[NSD];
708: PetscFunctionBeginUser;
709: PetscCall(DMGetCoordinatesLocal(dm_vp, &coor_l));
710: PetscCall(VecGetArrayRead(coor_l, &LA_coor));
712: PetscCall(DMGetLocalVector(dm_vp, &vp_l));
713: PetscCall(DMGlobalToLocalBegin(dm_vp, vp, INSERT_VALUES, vp_l));
714: PetscCall(DMGlobalToLocalEnd(dm_vp, vp, INSERT_VALUES, vp_l));
715: PetscCall(VecGetArrayRead(vp_l, &LA_vp));
717: PetscCall(DMDAGetElements(dm_vp, &nel, &npe, &element_list));
718: PetscCall(DMSwarmGetLocalSize(dm_mpoint, &npoints));
719: PetscCall(DMSwarmGetField(dm_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
720: PetscCall(DMSwarmGetField(dm_mpoint, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
721: for (p = 0; p < npoints; p++) {
722: PetscReal *coor_p;
723: PetscScalar vel_n[NSD * NODES_PER_EL], vel_p[NSD];
724: const PetscScalar *x0;
725: const PetscScalar *x2;
727: e = mpfield_cell[p];
728: coor_p = &mpfield_coor[NSD * p];
729: element = &element_list[NODES_PER_EL * e];
731: /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
732: x0 = &LA_coor[NSD * element[0]];
733: x2 = &LA_coor[NSD * element[2]];
735: dx[0] = x2[0] - x0[0];
736: dx[1] = x2[1] - x0[1];
738: xi_p[0] = 2.0 * (coor_p[0] - x0[0]) / dx[0] - 1.0;
739: xi_p[1] = 2.0 * (coor_p[1] - x0[1]) / dx[1] - 1.0;
740: 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);
741: 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);
742: 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);
743: 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);
745: /* evaluate basis functions */
746: EvaluateBasis_Q1(xi_p, Ni);
748: /* get cell nodal velocities */
749: for (i = 0; i < NODES_PER_EL; i++) {
750: PetscInt nid;
752: nid = element[i];
753: vel_n[NSD * i + 0] = LA_vp[(NSD + 1) * nid + 0];
754: vel_n[NSD * i + 1] = LA_vp[(NSD + 1) * nid + 1];
755: }
757: /* interpolate velocity */
758: vel_p[0] = vel_p[1] = 0.0;
759: for (i = 0; i < NODES_PER_EL; i++) {
760: vel_p[0] += Ni[i] * vel_n[NSD * i + 0];
761: vel_p[1] += Ni[i] * vel_n[NSD * i + 1];
762: }
764: coor_p[0] += dt * PetscRealPart(vel_p[0]);
765: coor_p[1] += dt * PetscRealPart(vel_p[1]);
766: }
768: PetscCall(DMSwarmRestoreField(dm_mpoint, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
769: PetscCall(DMSwarmRestoreField(dm_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
770: PetscCall(DMDARestoreElements(dm_vp, &nel, &npe, &element_list));
771: PetscCall(VecRestoreArrayRead(vp_l, &LA_vp));
772: PetscCall(DMRestoreLocalVector(dm_vp, &vp_l));
773: PetscCall(VecRestoreArrayRead(coor_l, &LA_coor));
774: PetscFunctionReturn(PETSC_SUCCESS);
775: }
777: PetscErrorCode MaterialPoint_Interpolate(DM dm, Vec eta_v, Vec rho_v, DM dm_quadrature)
778: {
779: Vec eta_l, rho_l;
780: PetscScalar *_eta_l, *_rho_l;
781: PetscInt nqp, npe, nel;
782: PetscScalar qp_xi[GAUSS_POINTS][NSD];
783: PetscScalar qp_weight[GAUSS_POINTS];
784: PetscInt q, k, e;
785: PetscScalar Ni[GAUSS_POINTS][NODES_PER_EL];
786: const PetscInt *element_list;
787: PetscReal *q_eta, *q_rhs;
789: PetscFunctionBeginUser;
790: /* define quadrature rule */
791: CreateGaussQuadrature(&nqp, qp_xi, qp_weight);
792: for (q = 0; q < nqp; q++) EvaluateBasis_Q1(qp_xi[q], Ni[q]);
794: PetscCall(DMGetLocalVector(dm, &eta_l));
795: PetscCall(DMGetLocalVector(dm, &rho_l));
797: PetscCall(DMGlobalToLocalBegin(dm, eta_v, INSERT_VALUES, eta_l));
798: PetscCall(DMGlobalToLocalEnd(dm, eta_v, INSERT_VALUES, eta_l));
799: PetscCall(DMGlobalToLocalBegin(dm, rho_v, INSERT_VALUES, rho_l));
800: PetscCall(DMGlobalToLocalEnd(dm, rho_v, INSERT_VALUES, rho_l));
802: PetscCall(VecGetArray(eta_l, &_eta_l));
803: PetscCall(VecGetArray(rho_l, &_rho_l));
805: PetscCall(DMSwarmGetField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
806: PetscCall(DMSwarmGetField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
808: PetscCall(DMDAGetElements(dm, &nel, &npe, &element_list));
809: for (e = 0; e < nel; e++) {
810: PetscScalar eta_field_e[NODES_PER_EL];
811: PetscScalar rho_field_e[NODES_PER_EL];
812: const PetscInt *element = &element_list[4 * e];
814: for (k = 0; k < NODES_PER_EL; k++) {
815: eta_field_e[k] = _eta_l[element[k]];
816: rho_field_e[k] = _rho_l[element[k]];
817: }
819: for (q = 0; q < nqp; q++) {
820: PetscScalar eta_q, rho_q;
822: eta_q = rho_q = 0.0;
823: for (k = 0; k < NODES_PER_EL; k++) {
824: eta_q += Ni[q][k] * eta_field_e[k];
825: rho_q += Ni[q][k] * rho_field_e[k];
826: }
828: q_eta[nqp * e + q] = PetscRealPart(eta_q);
829: q_rhs[nqp * e + q] = PetscRealPart(rho_q);
830: }
831: }
832: PetscCall(DMDARestoreElements(dm, &nel, &npe, &element_list));
834: PetscCall(DMSwarmRestoreField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
835: PetscCall(DMSwarmRestoreField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
837: PetscCall(VecRestoreArray(rho_l, &_rho_l));
838: PetscCall(VecRestoreArray(eta_l, &_eta_l));
839: PetscCall(DMRestoreLocalVector(dm, &rho_l));
840: PetscCall(DMRestoreLocalVector(dm, &eta_l));
841: PetscFunctionReturn(PETSC_SUCCESS);
842: }
844: static PetscErrorCode SolveTimeDepStokes(PetscInt mx, PetscInt my)
845: {
846: DM dm_stokes, dm_coeff;
847: PetscInt u_dof, p_dof, dof, stencil_width;
848: Mat A, B;
849: PetscInt nel_local;
850: Vec eta_v, rho_v;
851: Vec f, X;
852: KSP ksp;
853: PC pc;
854: char filename[PETSC_MAX_PATH_LEN];
855: DM dms_quadrature, dms_mpoint;
856: PetscInt nel, npe, npoints;
857: const PetscInt *element_list;
858: PetscInt tk, nt, dump_freq;
859: PetscReal dt, dt_max = 0.0;
860: PetscReal vx[2], vy[2], max_v = 0.0, max_v_step, dh;
861: const char *fieldnames[] = {"eta", "rho"};
862: Vec pfields[2];
863: PetscInt ppcell = 1;
864: PetscReal time, delta_eta = 1.0;
865: PetscBool randomize_coords = PETSC_FALSE;
866: PetscReal randomize_fac = 0.25;
867: PetscBool no_view = PETSC_FALSE;
868: PetscBool isbddc;
870: PetscFunctionBeginUser;
871: /*
872: Generate the DMDA for the velocity and pressure spaces.
873: We use Q1 elements for both fields.
874: The Q1 FE basis on a regular mesh has a 9-point stencil (DMDA_STENCIL_BOX)
875: The number of nodes in each direction is mx+1, my+1
876: */
877: u_dof = U_DOFS; /* Vx, Vy - velocities */
878: p_dof = P_DOFS; /* p - pressure */
879: dof = u_dof + p_dof;
880: stencil_width = 1;
881: 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));
882: PetscCall(DMDASetElementType(dm_stokes, DMDA_ELEMENT_Q1));
883: PetscCall(DMSetMatType(dm_stokes, MATAIJ));
884: PetscCall(DMSetFromOptions(dm_stokes));
885: PetscCall(DMSetUp(dm_stokes));
886: PetscCall(DMDASetFieldName(dm_stokes, 0, "ux"));
887: PetscCall(DMDASetFieldName(dm_stokes, 1, "uy"));
888: PetscCall(DMDASetFieldName(dm_stokes, 2, "p"));
890: /* unit box [0,0.9142] x [0,1] */
891: PetscCall(DMDASetUniformCoordinates(dm_stokes, 0.0, 0.9142, 0.0, 1.0, 0., 0.));
892: dh = 1.0 / ((PetscReal)(mx));
894: /* Get local number of elements */
895: {
896: PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));
898: nel_local = nel;
900: PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
901: }
903: /* Create DMDA for representing scalar fields */
904: PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 1, &dm_coeff));
906: /* Create the swarm for storing quadrature point values */
907: PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_quadrature));
908: PetscCall(DMSetType(dms_quadrature, DMSWARM));
909: PetscCall(DMSetDimension(dms_quadrature, 2));
910: PetscCall(PetscObjectSetName((PetscObject)dms_quadrature, "Quadrature Swarm"));
912: /* Register fields for viscosity and density on the quadrature points */
913: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "eta_q", 1, PETSC_REAL));
914: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "rho_q", 1, PETSC_REAL));
915: PetscCall(DMSwarmFinalizeFieldRegister(dms_quadrature));
916: PetscCall(DMSwarmSetLocalSizes(dms_quadrature, nel_local * GAUSS_POINTS, 0));
918: /* Create the material point swarm */
919: PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_mpoint));
920: PetscCall(DMSetType(dms_mpoint, DMSWARM));
921: PetscCall(DMSetDimension(dms_mpoint, 2));
922: PetscCall(PetscObjectSetName((PetscObject)dms_mpoint, "Material Point Swarm"));
924: /* Configure the material point swarm to be of type Particle-In-Cell */
925: PetscCall(DMSwarmSetType(dms_mpoint, DMSWARM_PIC));
927: /*
928: Specify the DM to use for point location and projections
929: within the context of a PIC scheme
930: */
931: PetscCall(DMSwarmSetCellDM(dms_mpoint, dm_coeff));
933: /* Register fields for viscosity and density */
934: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "eta", 1, PETSC_REAL));
935: PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "rho", 1, PETSC_REAL));
936: PetscCall(DMSwarmFinalizeFieldRegister(dms_mpoint));
938: PetscCall(PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL));
939: PetscCall(DMSwarmSetLocalSizes(dms_mpoint, nel_local * ppcell, 100));
941: /*
942: Layout the material points in space using the cell DM.
943: Particle coordinates are defined by cell wise using different methods.
944: - DMSWARMPIC_LAYOUT_GAUSS defines particles coordinates at the positions
945: corresponding to a Gauss quadrature rule with
946: ppcell points in each direction.
947: - DMSWARMPIC_LAYOUT_REGULAR defines particle coordinates at the centoid of
948: ppcell x ppcell quadralaterals defined within the
949: reference element.
950: - DMSWARMPIC_LAYOUT_SUBDIVISION defines particles coordinates at the centroid
951: of each quadralateral obtained by sub-dividing
952: the reference element cell ppcell times.
953: */
954: PetscCall(DMSwarmInsertPointsUsingCellDM(dms_mpoint, DMSWARMPIC_LAYOUT_SUBDIVISION, ppcell));
956: /*
957: Defne a high resolution layer of material points across the material interface
958: */
959: {
960: PetscInt npoints_dir_x[2];
961: PetscReal min[2], max[2];
963: npoints_dir_x[0] = (PetscInt)(0.9142 / (0.05 * dh));
964: npoints_dir_x[1] = (PetscInt)((0.25 - 0.15) / (0.05 * dh));
965: min[0] = 0.0;
966: max[0] = 0.9142;
967: min[1] = 0.05;
968: max[1] = 0.35;
969: PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
970: }
972: /*
973: Define a high resolution layer of material points near the surface of the domain
974: to deal with weakly compressible Q1-Q1 elements. These elements "self compact"
975: when applied to buouyancy driven flow. The error in div(u) is O(h).
976: */
977: {
978: PetscInt npoints_dir_x[2];
979: PetscReal min[2], max[2];
981: npoints_dir_x[0] = (PetscInt)(0.9142 / (0.25 * dh));
982: npoints_dir_x[1] = (PetscInt)(3.0 * dh / (0.25 * dh));
983: min[0] = 0.0;
984: max[0] = 0.9142;
985: min[1] = 1.0 - 3.0 * dh;
986: max[1] = 1.0 - 0.0001;
987: PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
988: }
990: PetscCall(DMView(dms_mpoint, PETSC_VIEWER_STDOUT_WORLD));
992: /* Define initial material properties on each particle in the material point swarm */
993: PetscCall(PetscOptionsGetReal(NULL, NULL, "-delta_eta", &delta_eta, NULL));
994: PetscCall(PetscOptionsGetBool(NULL, NULL, "-randomize_coords", &randomize_coords, NULL));
995: PetscCall(PetscOptionsGetReal(NULL, NULL, "-randomize_fac", &randomize_fac, NULL));
996: PetscCheck(randomize_fac <= 1.0, PETSC_COMM_WORLD, PETSC_ERR_USER, "The value of -randomize_fac should be <= 1.0");
997: {
998: PetscReal *array_x, *array_e, *array_r;
999: PetscInt p;
1000: PetscRandom r;
1001: PetscMPIInt rank;
1003: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
1005: PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
1006: PetscCall(PetscRandomSetInterval(r, -randomize_fac * dh, randomize_fac * dh));
1007: PetscCall(PetscRandomSetSeed(r, (unsigned long)rank));
1008: PetscCall(PetscRandomSeed(r));
1010: PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));
1012: /*
1013: Fetch the registered data from the material point DMSwarm.
1014: The fields "eta" and "rho" were registered by this example.
1015: The field identified by the variable DMSwarmPICField_coor
1016: was registered by the DMSwarm implementation when the function
1017: DMSwarmSetType(dms_mpoint,DMSWARM_PIC)
1018: was called. The returned array defines the coordinates of each
1019: material point in the point swarm.
1020: */
1021: PetscCall(DMSwarmGetField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));
1022: PetscCall(DMSwarmGetField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
1023: PetscCall(DMSwarmGetField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));
1025: PetscCall(DMSwarmGetLocalSize(dms_mpoint, &npoints));
1026: for (p = 0; p < npoints; p++) {
1027: PetscReal x_p[2], rr[2];
1029: if (randomize_coords) {
1030: PetscCall(PetscRandomGetValueReal(r, &rr[0]));
1031: PetscCall(PetscRandomGetValueReal(r, &rr[1]));
1032: array_x[2 * p + 0] += rr[0];
1033: array_x[2 * p + 1] += rr[1];
1034: }
1036: /* Get the coordinates of point, p */
1037: x_p[0] = array_x[2 * p + 0];
1038: x_p[1] = array_x[2 * p + 1];
1040: if (x_p[1] < (0.2 + 0.02 * PetscCosReal(PETSC_PI * x_p[0] / 0.9142))) {
1041: /* Material properties below the interface */
1042: array_e[p] = 1.0 * (1.0 / delta_eta);
1043: array_r[p] = 0.0;
1044: } else {
1045: /* Material properties above the interface */
1046: array_e[p] = 1.0;
1047: array_r[p] = 1.0;
1048: }
1049: }
1051: /*
1052: Restore the fetched data fields from the material point DMSwarm.
1053: Calling the Restore function invalidates the points array_r, array_e, array_x
1054: by setting them to NULL.
1055: */
1056: PetscCall(DMSwarmRestoreField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));
1057: PetscCall(DMSwarmRestoreField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
1058: PetscCall(DMSwarmRestoreField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));
1060: PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
1061: PetscCall(PetscRandomDestroy(&r));
1062: }
1064: /*
1065: If the particle coordinates where randomly shifted, they may have crossed into another
1066: element, or into another sub-domain. To account for this we call the Migrate function.
1067: */
1068: if (randomize_coords) PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));
1070: PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL));
1071: if (!no_view) PetscCall(DMSwarmViewXDMF(dms_mpoint, "ic_coeff_dms.xmf"));
1073: /* project the swarm properties */
1074: PetscCall(DMCreateGlobalVector(dm_coeff, &pfields[0]));
1075: PetscCall(DMCreateGlobalVector(dm_coeff, &pfields[1]));
1076: PetscCall(DMSwarmProjectFields(dms_mpoint, NULL, 2, fieldnames, pfields, SCATTER_FORWARD));
1077: eta_v = pfields[0];
1078: rho_v = pfields[1];
1079: PetscCall(PetscObjectSetName((PetscObject)eta_v, "eta"));
1080: PetscCall(PetscObjectSetName((PetscObject)rho_v, "rho"));
1081: PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));
1083: /* view projected coefficients eta and rho */
1084: if (!no_view) {
1085: PetscViewer viewer;
1087: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1088: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1089: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1090: PetscCall(PetscViewerFileSetName(viewer, "ic_coeff_dmda.vts"));
1091: PetscCall(VecView(eta_v, viewer));
1092: PetscCall(VecView(rho_v, viewer));
1093: PetscCall(PetscViewerDestroy(&viewer));
1094: }
1096: PetscCall(DMCreateMatrix(dm_stokes, &A));
1097: PetscCall(DMCreateMatrix(dm_stokes, &B));
1098: PetscCall(DMCreateGlobalVector(dm_stokes, &f));
1099: PetscCall(DMCreateGlobalVector(dm_stokes, &X));
1101: PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
1102: PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
1103: PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));
1105: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
1106: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));
1108: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
1109: PetscCall(KSPSetOptionsPrefix(ksp, "stokes_"));
1110: PetscCall(KSPSetDM(ksp, dm_stokes));
1111: PetscCall(KSPSetDMActive(ksp, PETSC_FALSE));
1112: PetscCall(KSPSetOperators(ksp, A, B));
1113: PetscCall(KSPSetFromOptions(ksp));
1114: PetscCall(KSPGetPC(ksp, &pc));
1115: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &isbddc));
1116: if (isbddc) PetscCall(KSPSetOperators(ksp, A, A));
1118: /* Define u-v-p indices for fieldsplit */
1119: {
1120: PC pc;
1121: const PetscInt ufields[] = {0, 1}, pfields[1] = {2};
1123: PetscCall(KSPGetPC(ksp, &pc));
1124: PetscCall(PCFieldSplitSetBlockSize(pc, 3));
1125: PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
1126: PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
1127: }
1129: /* If using a fieldsplit preconditioner, attach a DMDA to the velocity split so that geometric multigrid can be used */
1130: {
1131: PC pc, pc_u;
1132: KSP *sub_ksp, ksp_u;
1133: PetscInt nsplits;
1134: DM dm_u;
1135: PetscBool is_pcfs;
1137: PetscCall(KSPGetPC(ksp, &pc));
1139: is_pcfs = PETSC_FALSE;
1140: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs));
1142: if (is_pcfs) {
1143: PetscCall(KSPSetUp(ksp));
1144: PetscCall(KSPGetPC(ksp, &pc));
1145: PetscCall(PCFieldSplitGetSubKSP(pc, &nsplits, &sub_ksp));
1146: ksp_u = sub_ksp[0];
1147: PetscCall(PetscFree(sub_ksp));
1149: if (nsplits == 2) {
1150: PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 2, &dm_u));
1152: PetscCall(KSPSetDM(ksp_u, dm_u));
1153: PetscCall(KSPSetDMActive(ksp_u, PETSC_FALSE));
1154: PetscCall(DMDestroy(&dm_u));
1156: /* enforce galerkin coarse grids be used */
1157: PetscCall(KSPGetPC(ksp_u, &pc_u));
1158: PetscCall(PCMGSetGalerkin(pc_u, PC_MG_GALERKIN_PMAT));
1159: }
1160: }
1161: }
1163: dump_freq = 10;
1164: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dump_freq", &dump_freq, NULL));
1165: nt = 10;
1166: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nt", &nt, NULL));
1167: time = 0.0;
1168: for (tk = 1; tk <= nt; tk++) {
1169: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... assemble\n"));
1170: PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
1171: PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
1172: PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));
1174: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... bc imposition\n"));
1175: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
1176: PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));
1178: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... solve\n"));
1179: PetscCall(KSPSetOperators(ksp, A, isbddc ? A : B));
1180: PetscCall(KSPSolve(ksp, f, X));
1182: PetscCall(VecStrideMax(X, 0, NULL, &vx[1]));
1183: PetscCall(VecStrideMax(X, 1, NULL, &vy[1]));
1184: PetscCall(VecStrideMin(X, 0, NULL, &vx[0]));
1185: PetscCall(VecStrideMin(X, 1, NULL, &vy[0]));
1187: max_v_step = PetscMax(vx[0], vx[1]);
1188: max_v_step = PetscMax(max_v_step, vy[0]);
1189: max_v_step = PetscMax(max_v_step, vy[1]);
1190: max_v = PetscMax(max_v, max_v_step);
1192: dt_max = 2.0;
1193: dt = 0.5 * (dh / max_v_step);
1194: 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));
1195: dt = PetscMin(dt_max, dt);
1197: /* advect */
1198: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... advect\n"));
1199: PetscCall(MaterialPoint_AdvectRK1(dm_stokes, X, dt, dms_mpoint));
1201: /* migrate */
1202: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... migrate\n"));
1203: PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));
1205: /* update cell population */
1206: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... populate cells\n"));
1207: PetscCall(MaterialPoint_PopulateCell(dm_stokes, dms_mpoint));
1209: /* update coefficients on quadrature points */
1210: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... project\n"));
1211: PetscCall(DMSwarmProjectFields(dms_mpoint, NULL, 2, fieldnames, pfields, SCATTER_FORWARD));
1212: eta_v = pfields[0];
1213: rho_v = pfields[1];
1214: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... interp\n"));
1215: PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));
1217: if (tk % dump_freq == 0) {
1218: PetscViewer viewer;
1220: PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... write XDMF, VTS\n"));
1221: PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_coeff_dms.xmf", tk));
1222: PetscCall(DMSwarmViewXDMF(dms_mpoint, filename));
1224: PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_vp_dm.vts", tk));
1225: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1226: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1227: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1228: PetscCall(PetscViewerFileSetName(viewer, filename));
1229: PetscCall(VecView(X, viewer));
1230: PetscCall(PetscViewerDestroy(&viewer));
1231: }
1232: time += dt;
1233: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step %" PetscInt_FMT " : time %1.2e \n", tk, (double)time));
1234: }
1236: PetscCall(KSPDestroy(&ksp));
1237: PetscCall(VecDestroy(&X));
1238: PetscCall(VecDestroy(&f));
1239: PetscCall(MatDestroy(&A));
1240: PetscCall(MatDestroy(&B));
1241: PetscCall(VecDestroy(&eta_v));
1242: PetscCall(VecDestroy(&rho_v));
1244: PetscCall(DMDestroy(&dms_mpoint));
1245: PetscCall(DMDestroy(&dms_quadrature));
1246: PetscCall(DMDestroy(&dm_coeff));
1247: PetscCall(DMDestroy(&dm_stokes));
1248: PetscFunctionReturn(PETSC_SUCCESS);
1249: }
1251: /*
1252: <sequential run>
1253: ./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
1254: */
1255: int main(int argc, char **args)
1256: {
1257: PetscInt mx, my;
1258: PetscBool set = PETSC_FALSE;
1260: PetscFunctionBeginUser;
1261: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
1262: mx = my = 10;
1263: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1264: PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1265: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mxy", &mx, &set));
1266: if (set) my = mx;
1267: PetscCall(SolveTimeDepStokes(mx, my));
1268: PetscCall(PetscFinalize());
1269: return 0;
1270: }
1272: /* -------------------------- helpers for boundary conditions -------------------------------- */
1273: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
1274: {
1275: DM cda;
1276: Vec coords;
1277: PetscInt si, sj, nx, ny, i, j;
1278: PetscInt M, N;
1279: DMDACoor2d **_coords;
1280: const PetscInt *g_idx;
1281: PetscInt *bc_global_ids;
1282: PetscScalar *bc_vals;
1283: PetscInt nbcs;
1284: PetscInt n_dofs;
1285: ISLocalToGlobalMapping ltogm;
1287: PetscFunctionBeginUser;
1288: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1289: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1291: PetscCall(DMGetCoordinateDM(da, &cda));
1292: PetscCall(DMGetCoordinatesLocal(da, &coords));
1293: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1294: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1295: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1297: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1298: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1300: /* init the entries to -1 so VecSetValues will ignore them */
1301: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1303: i = nx - 1;
1304: for (j = 0; j < ny; j++) {
1305: PetscInt local_id;
1307: local_id = i + j * nx;
1309: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1311: bc_vals[j] = 0.0;
1312: }
1313: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1314: nbcs = 0;
1315: if ((si + nx) == (M)) nbcs = ny;
1317: if (b) {
1318: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1319: PetscCall(VecAssemblyBegin(b));
1320: PetscCall(VecAssemblyEnd(b));
1321: }
1322: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1324: PetscCall(PetscFree(bc_vals));
1325: PetscCall(PetscFree(bc_global_ids));
1327: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1328: PetscFunctionReturn(PETSC_SUCCESS);
1329: }
1331: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
1332: {
1333: DM cda;
1334: Vec coords;
1335: PetscInt si, sj, nx, ny, i, j;
1336: PetscInt M, N;
1337: DMDACoor2d **_coords;
1338: const PetscInt *g_idx;
1339: PetscInt *bc_global_ids;
1340: PetscScalar *bc_vals;
1341: PetscInt nbcs;
1342: PetscInt n_dofs;
1343: ISLocalToGlobalMapping ltogm;
1345: PetscFunctionBeginUser;
1346: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1347: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1349: PetscCall(DMGetCoordinateDM(da, &cda));
1350: PetscCall(DMGetCoordinatesLocal(da, &coords));
1351: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1352: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1353: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1355: PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1356: PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));
1358: /* init the entries to -1 so VecSetValues will ignore them */
1359: for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;
1361: i = 0;
1362: for (j = 0; j < ny; j++) {
1363: PetscInt local_id;
1365: local_id = i + j * nx;
1367: bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];
1369: bc_vals[j] = 0.0;
1370: }
1371: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1372: nbcs = 0;
1373: if (si == 0) nbcs = ny;
1375: if (b) {
1376: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1377: PetscCall(VecAssemblyBegin(b));
1378: PetscCall(VecAssemblyEnd(b));
1379: }
1381: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1383: PetscCall(PetscFree(bc_vals));
1384: PetscCall(PetscFree(bc_global_ids));
1386: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1387: PetscFunctionReturn(PETSC_SUCCESS);
1388: }
1390: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
1391: {
1392: DM cda;
1393: Vec coords;
1394: PetscInt si, sj, nx, ny, i, j;
1395: PetscInt M, N;
1396: DMDACoor2d **_coords;
1397: const PetscInt *g_idx;
1398: PetscInt *bc_global_ids;
1399: PetscScalar *bc_vals;
1400: PetscInt nbcs;
1401: PetscInt n_dofs;
1402: ISLocalToGlobalMapping ltogm;
1404: PetscFunctionBeginUser;
1405: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1406: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1408: PetscCall(DMGetCoordinateDM(da, &cda));
1409: PetscCall(DMGetCoordinatesLocal(da, &coords));
1410: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1411: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1412: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1414: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1415: PetscCall(PetscMalloc1(nx, &bc_vals));
1417: /* init the entries to -1 so VecSetValues will ignore them */
1418: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1420: j = ny - 1;
1421: for (i = 0; i < nx; i++) {
1422: PetscInt local_id;
1424: local_id = i + j * nx;
1426: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1428: bc_vals[i] = 0.0;
1429: }
1430: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1431: nbcs = 0;
1432: if ((sj + ny) == (N)) nbcs = nx;
1434: if (b) {
1435: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1436: PetscCall(VecAssemblyBegin(b));
1437: PetscCall(VecAssemblyEnd(b));
1438: }
1439: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, NULL, NULL));
1441: PetscCall(PetscFree(bc_vals));
1442: PetscCall(PetscFree(bc_global_ids));
1444: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1445: PetscFunctionReturn(PETSC_SUCCESS);
1446: }
1448: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
1449: {
1450: DM cda;
1451: Vec coords;
1452: PetscInt si, sj, nx, ny, i, j;
1453: PetscInt M, N;
1454: DMDACoor2d **_coords;
1455: const PetscInt *g_idx;
1456: PetscInt *bc_global_ids;
1457: PetscScalar *bc_vals;
1458: PetscInt nbcs;
1459: PetscInt n_dofs;
1460: ISLocalToGlobalMapping ltogm;
1462: PetscFunctionBeginUser;
1463: PetscCall(DMGetLocalToGlobalMapping(da, <ogm));
1464: PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));
1466: PetscCall(DMGetCoordinateDM(da, &cda));
1467: PetscCall(DMGetCoordinatesLocal(da, &coords));
1468: PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1469: PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1470: PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));
1472: PetscCall(PetscMalloc1(nx, &bc_global_ids));
1473: PetscCall(PetscMalloc1(nx, &bc_vals));
1475: /* init the entries to -1 so VecSetValues will ignore them */
1476: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1478: j = 0;
1479: for (i = 0; i < nx; i++) {
1480: PetscInt local_id;
1482: local_id = i + j * nx;
1484: bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];
1486: bc_vals[i] = 0.0;
1487: }
1488: PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1489: nbcs = 0;
1490: if (sj == 0) nbcs = nx;
1492: if (b) {
1493: PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1494: PetscCall(VecAssemblyBegin(b));
1495: PetscCall(VecAssemblyEnd(b));
1496: }
1497: if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));
1499: PetscCall(PetscFree(bc_vals));
1500: PetscCall(PetscFree(bc_global_ids));
1502: PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1503: PetscFunctionReturn(PETSC_SUCCESS);
1504: }
1506: /*
1507: Impose free slip boundary conditions on the left/right faces: u_i n_i = 0, tau_{ij} t_j = 0
1508: Impose no slip boundray conditions on the top/bottom faces: u_i n_i = 0, u_i t_i = 0
1509: */
1510: static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes, Mat A, Vec f)
1511: {
1512: PetscFunctionBeginUser;
1513: PetscCall(BCApplyZero_NORTH(dm_stokes, 0, A, f));
1514: PetscCall(BCApplyZero_NORTH(dm_stokes, 1, A, f));
1515: PetscCall(BCApplyZero_EAST(dm_stokes, 0, A, f));
1516: PetscCall(BCApplyZero_SOUTH(dm_stokes, 0, A, f));
1517: PetscCall(BCApplyZero_SOUTH(dm_stokes, 1, A, f));
1518: PetscCall(BCApplyZero_WEST(dm_stokes, 0, A, f));
1519: PetscFunctionReturn(PETSC_SUCCESS);
1520: }
1522: /*TEST
1524: test:
1525: suffix: 1
1526: args: -no_view
1527: requires: !complex double
1528: filter: grep -v atomic
1529: filter_output: grep -v atomic
1530: test:
1531: suffix: 1_matis
1532: requires: !complex double
1533: args: -no_view -dm_mat_type is
1534: filter: grep -v atomic
1535: filter_output: grep -v atomic
1536: testset:
1537: nsize: 4
1538: requires: !complex double
1539: 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
1540: filter: grep -v atomic
1541: filter_output: grep -v atomic
1542: test:
1543: suffix: fetidp
1544: args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
1545: test:
1546: suffix: fetidp_lumped
1547: 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
1548: test:
1549: suffix: fetidp_saddlepoint
1550: 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
1551: test:
1552: suffix: fetidp_saddlepoint_lumped
1553: args: -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_norm_type natural -stokes_fetidp_pc_fieldsplit_schur_fact_type diag -stokes_fetidp_fieldsplit_p_pc_type bjacobi -stokes_fetidp_fieldsplit_lag_ksp_type preonly -stokes_fetidp_fieldsplit_p_ksp_type preonly -stokes_ksp_fetidp_pressure_field 2 -stokes_fetidp_pc_fieldsplit_schur_scale -1 -stokes_fetidp_bddc_pc_bddc_dirichlet_pc_type none -stokes_fetidp_bddc_pc_bddc_switch_static -stokes_fetidp_pc_lumped
1554: TEST*/