Actual source code: ex70.c

  1: static char help[] = "------------------------------------------------------------------------------------------------------------------------------ \n\
  2:   Solves the time-dependent incompressible, variable viscosity Stokes equation in 2D driven by buouyancy variations. \n\
  3:   Time-dependence is introduced by evolving the density (rho) and viscosity (eta) according to \n\
  4:     D \\rho / Dt = 0    and    D \\eta / Dt = 0 \n\
  5:   The Stokes problem is discretized using Q1-Q1 finite elements, stabilized with Bochev's polynomial projection method. \n\
  6:   The hyperbolic evolution equation for density is discretized using a variant of the Particle-In-Cell (PIC) method. \n\
  7:   The DMDA object is used to define the FE problem, whilst DMSwarm provides support for the PIC method. \n\
  8:   Material points (particles) store density and viscosity. The particles are advected with the fluid velocity using RK1. \n\
  9:   At each time step, the value of density and viscosity stored on each particle are projected into a Q1 function space \n\
 10:   and then interpolated onto the Gauss quadrature points. \n\
 11:   The model problem defined in this example is the iso-viscous Rayleigh-Taylor instability (case 1a) from: \n\
 12:     \"A comparison of methods for the modeling of thermochemical convection\" \n\
 13:     P.E. van Keken, S.D. King, H. Schmeling, U.R. Christensen, D. Neumeister and M.-P. Doin, \n\
 14:     Journal of Geophysical Research, vol 102 (B10), 477--499 (1997) \n\
 15:   Note that whilst the model problem defined is for an iso-viscous, the implementation in this example supports \n\
 16:   variable viscoity formulations. \n\
 17:   This example is based on src/ksp/ksp/tutorials/ex43.c \n\
 18:   Options: \n\
 19:     -mx        : Number of elements in the x-direction \n\
 20:     -my        : Number of elements in the y-direction \n\
 21:     -mxy       : Number of elements in the x- and y-directions \n\
 22:     -nt        : Number of time steps \n\
 23:     -dump_freq : Frequency of output file creation \n\
 24:     -ppcell    : Number of times the reference cell is sub-divided \n\
 25:     -randomize_coords : Apply a random shift to each particle coordinate in the range [-fac*dh,0.fac*dh] \n\
 26:     -randomize_fac    : Set the scaling factor for the random shift (default = 0.25)\n";

 28: /* Contributed by Dave May */

 30: #include <petscksp.h>
 31: #include <petscdm.h>
 32: #include <petscdmda.h>
 33: #include <petscdmswarm.h>
 34: #include <petsc/private/dmimpl.h>

 36: static PetscErrorCode DMDAApplyBoundaryConditions(DM, Mat, Vec);

 38: #define NSD          2 /* number of spatial dimensions */
 39: #define NODES_PER_EL 4 /* nodes per element */
 40: #define U_DOFS       2 /* degrees of freedom per velocity node */
 41: #define P_DOFS       1 /* degrees of freedom per pressure node */
 42: #define GAUSS_POINTS 4

 44: static void EvaluateBasis_Q1(PetscScalar _xi[], PetscScalar N[])
 45: {
 46:   PetscScalar xi  = _xi[0];
 47:   PetscScalar eta = _xi[1];

 49:   N[0] = 0.25 * (1.0 - xi) * (1.0 - eta);
 50:   N[1] = 0.25 * (1.0 + xi) * (1.0 - eta);
 51:   N[2] = 0.25 * (1.0 + xi) * (1.0 + eta);
 52:   N[3] = 0.25 * (1.0 - xi) * (1.0 + eta);
 53: }

 55: static void EvaluateBasisDerivatives_Q1(PetscScalar _xi[], PetscScalar dN[][NODES_PER_EL])
 56: {
 57:   PetscScalar xi  = _xi[0];
 58:   PetscScalar eta = _xi[1];

 60:   dN[0][0] = -0.25 * (1.0 - eta);
 61:   dN[0][1] = 0.25 * (1.0 - eta);
 62:   dN[0][2] = 0.25 * (1.0 + eta);
 63:   dN[0][3] = -0.25 * (1.0 + eta);

 65:   dN[1][0] = -0.25 * (1.0 - xi);
 66:   dN[1][1] = -0.25 * (1.0 + xi);
 67:   dN[1][2] = 0.25 * (1.0 + xi);
 68:   dN[1][3] = 0.25 * (1.0 - xi);
 69: }

 71: static void EvaluateDerivatives(PetscScalar dN[][NODES_PER_EL], PetscScalar dNx[][NODES_PER_EL], PetscScalar coords[], PetscScalar *det_J)
 72: {
 73:   PetscScalar J00, J01, J10, J11, J;
 74:   PetscScalar iJ00, iJ01, iJ10, iJ11;
 75:   PetscInt    i;

 77:   J00 = J01 = J10 = J11 = 0.0;
 78:   for (i = 0; i < NODES_PER_EL; i++) {
 79:     PetscScalar cx = coords[2 * i];
 80:     PetscScalar cy = coords[2 * i + 1];

 82:     J00 += dN[0][i] * cx; /* J_xx = dx/dxi */
 83:     J01 += dN[0][i] * cy; /* J_xy = dy/dxi */
 84:     J10 += dN[1][i] * cx; /* J_yx = dx/deta */
 85:     J11 += dN[1][i] * cy; /* J_yy = dy/deta */
 86:   }
 87:   J = (J00 * J11) - (J01 * J10);

 89:   iJ00 = J11 / J;
 90:   iJ01 = -J01 / J;
 91:   iJ10 = -J10 / J;
 92:   iJ11 = J00 / J;

 94:   for (i = 0; i < NODES_PER_EL; i++) {
 95:     dNx[0][i] = dN[0][i] * iJ00 + dN[1][i] * iJ01;
 96:     dNx[1][i] = dN[0][i] * iJ10 + dN[1][i] * iJ11;
 97:   }

 99:   if (det_J) *det_J = J;
100: }

102: static void CreateGaussQuadrature(PetscInt *ngp, PetscScalar gp_xi[][2], PetscScalar gp_weight[])
103: {
104:   *ngp         = 4;
105:   gp_xi[0][0]  = -0.57735026919;
106:   gp_xi[0][1]  = -0.57735026919;
107:   gp_xi[1][0]  = -0.57735026919;
108:   gp_xi[1][1]  = 0.57735026919;
109:   gp_xi[2][0]  = 0.57735026919;
110:   gp_xi[2][1]  = 0.57735026919;
111:   gp_xi[3][0]  = 0.57735026919;
112:   gp_xi[3][1]  = -0.57735026919;
113:   gp_weight[0] = 1.0;
114:   gp_weight[1] = 1.0;
115:   gp_weight[2] = 1.0;
116:   gp_weight[3] = 1.0;
117: }

119: static PetscErrorCode DMDAGetElementEqnums_up(const PetscInt element[], PetscInt s_u[], PetscInt s_p[])
120: {
121:   PetscInt i;
122:   PetscFunctionBeginUser;
123:   for (i = 0; i < NODES_PER_EL; i++) {
124:     /* velocity */
125:     s_u[NSD * i + 0] = 3 * element[i];
126:     s_u[NSD * i + 1] = 3 * element[i] + 1;
127:     /* pressure */
128:     s_p[i] = 3 * element[i] + 2;
129:   }
130:   PetscFunctionReturn(PETSC_SUCCESS);
131: }

133: static PetscInt map_wIwDI_uJuDJ(PetscInt wi, PetscInt wd, PetscInt w_NPE, PetscInt w_dof, PetscInt ui, PetscInt ud, PetscInt u_NPE, PetscInt u_dof)
134: {
135:   PetscInt ij, r, c, nc;

137:   nc = u_NPE * u_dof;
138:   r  = w_dof * wi + wd;
139:   c  = u_dof * ui + ud;
140:   ij = r * nc + c;
141:   return (ij);
142: }

144: static void BForm_DivT(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
145: {
146:   PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
147:   PetscScalar GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
148:   PetscScalar J_p, tildeD[3];
149:   PetscScalar B[3][U_DOFS * NODES_PER_EL];
150:   PetscInt    p, i, j, k, ngp;

152:   /* define quadrature rule */
153:   CreateGaussQuadrature(&ngp, gp_xi, gp_weight);

155:   /* evaluate bilinear form */
156:   for (p = 0; p < ngp; p++) {
157:     EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
158:     EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);

160:     for (i = 0; i < NODES_PER_EL; i++) {
161:       PetscScalar d_dx_i = GNx_p[0][i];
162:       PetscScalar d_dy_i = GNx_p[1][i];

164:       B[0][2 * i]     = d_dx_i;
165:       B[0][2 * i + 1] = 0.0;
166:       B[1][2 * i]     = 0.0;
167:       B[1][2 * i + 1] = d_dy_i;
168:       B[2][2 * i]     = d_dy_i;
169:       B[2][2 * i + 1] = d_dx_i;
170:     }

172:     tildeD[0] = 2.0 * gp_weight[p] * J_p * eta[p];
173:     tildeD[1] = 2.0 * gp_weight[p] * J_p * eta[p];
174:     tildeD[2] = gp_weight[p] * J_p * eta[p];

176:     /* form Bt tildeD B */
177:     /*
178:     Ke_ij = Bt_ik . D_kl . B_lj
179:           = B_ki . D_kl . B_lj
180:           = B_ki . D_kk . B_kj
181:     */
182:     for (i = 0; i < 8; i++) {
183:       for (j = 0; j < 8; j++) {
184:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
185:           Ke[i + 8 * j] += B[k][i] * tildeD[k] * B[k][j];
186:         }
187:       }
188:     }
189:   }
190: }

192: static void BForm_Grad(PetscScalar Ke[], PetscScalar coords[])
193: {
194:   PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
195:   PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
196:   PetscScalar J_p, fac;
197:   PetscInt    p, i, j, di, ngp;

199:   /* define quadrature rule */
200:   CreateGaussQuadrature(&ngp, gp_xi, gp_weight);

202:   /* evaluate bilinear form */
203:   for (p = 0; p < ngp; p++) {
204:     EvaluateBasis_Q1(gp_xi[p], Ni_p);
205:     EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
206:     EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
207:     fac = gp_weight[p] * J_p;

209:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
210:       for (di = 0; di < NSD; di++) {     /* u dofs */
211:         for (j = 0; j < 4; j++) {        /* p nodes, p dofs = 1 (ie no loop) */
212:           PetscInt IJ;
213:           IJ = map_wIwDI_uJuDJ(i, di, NODES_PER_EL, 2, j, 0, NODES_PER_EL, 1);

215:           Ke[IJ] -= GNx_p[di][i] * Ni_p[j] * fac;
216:         }
217:       }
218:     }
219:   }
220: }

222: static void BForm_Div(PetscScalar De[], PetscScalar coords[])
223: {
224:   PetscScalar Ge[U_DOFS * NODES_PER_EL * P_DOFS * NODES_PER_EL];
225:   PetscInt    i, j, nr_g, nc_g;

227:   PetscCallAbort(PETSC_COMM_SELF, PetscMemzero(Ge, sizeof(Ge)));
228:   BForm_Grad(Ge, coords);

230:   nr_g = U_DOFS * NODES_PER_EL;
231:   nc_g = P_DOFS * NODES_PER_EL;

233:   for (i = 0; i < nr_g; i++) {
234:     for (j = 0; j < nc_g; j++) De[nr_g * j + i] = Ge[nc_g * i + j];
235:   }
236: }

238: static void BForm_Stabilisation(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
239: {
240:   PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
241:   PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
242:   PetscScalar J_p, fac, eta_avg;
243:   PetscInt    p, i, j, ngp;

245:   /* define quadrature rule */
246:   CreateGaussQuadrature(&ngp, gp_xi, gp_weight);

248:   /* evaluate bilinear form */
249:   for (p = 0; p < ngp; p++) {
250:     EvaluateBasis_Q1(gp_xi[p], Ni_p);
251:     EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
252:     EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
253:     fac = gp_weight[p] * J_p;

255:     for (i = 0; i < NODES_PER_EL; i++) {
256:       for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * (Ni_p[i] * Ni_p[j] - 0.0625);
257:     }
258:   }

260:   /* scale */
261:   eta_avg = 0.0;
262:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
263:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
264:   fac     = 1.0 / eta_avg;
265:   for (i = 0; i < NODES_PER_EL; i++) {
266:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] = fac * Ke[NODES_PER_EL * i + j];
267:   }
268: }

270: static void BForm_ScaledMassMatrix(PetscScalar Ke[], PetscScalar coords[], PetscScalar eta[])
271: {
272:   PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
273:   PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
274:   PetscScalar J_p, fac, eta_avg;
275:   PetscInt    p, i, j, ngp;

277:   /* define quadrature rule */
278:   CreateGaussQuadrature(&ngp, gp_xi, gp_weight);

280:   /* evaluate bilinear form */
281:   for (p = 0; p < ngp; p++) {
282:     EvaluateBasis_Q1(gp_xi[p], Ni_p);
283:     EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
284:     EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
285:     fac = gp_weight[p] * J_p;

287:     for (i = 0; i < NODES_PER_EL; i++) {
288:       for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] -= fac * Ni_p[i] * Ni_p[j];
289:     }
290:   }

292:   /* scale */
293:   eta_avg = 0.0;
294:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
295:   eta_avg = (1.0 / ((PetscScalar)ngp)) * eta_avg;
296:   fac     = 1.0 / eta_avg;
297:   for (i = 0; i < NODES_PER_EL; i++) {
298:     for (j = 0; j < NODES_PER_EL; j++) Ke[NODES_PER_EL * i + j] *= fac;
299:   }
300: }

302: static void LForm_MomentumRHS(PetscScalar Fe[], PetscScalar coords[], PetscScalar fx[], PetscScalar fy[])
303: {
304:   PetscScalar gp_xi[GAUSS_POINTS][NSD], gp_weight[GAUSS_POINTS];
305:   PetscScalar Ni_p[NODES_PER_EL], GNi_p[NSD][NODES_PER_EL], GNx_p[NSD][NODES_PER_EL];
306:   PetscScalar J_p, fac;
307:   PetscInt    p, i, ngp;

309:   /* define quadrature rule */
310:   CreateGaussQuadrature(&ngp, gp_xi, gp_weight);

312:   /* evaluate linear form */
313:   for (p = 0; p < ngp; p++) {
314:     EvaluateBasis_Q1(gp_xi[p], Ni_p);
315:     EvaluateBasisDerivatives_Q1(gp_xi[p], GNi_p);
316:     EvaluateDerivatives(GNi_p, GNx_p, coords, &J_p);
317:     fac = gp_weight[p] * J_p;

319:     for (i = 0; i < NODES_PER_EL; i++) {
320:       Fe[NSD * i] = 0.0;
321:       Fe[NSD * i + 1] -= fac * Ni_p[i] * fy[p];
322:     }
323:   }
324: }

326: static PetscErrorCode GetElementCoords(const PetscScalar _coords[], const PetscInt e2n[], PetscScalar el_coords[])
327: {
328:   PetscInt i, d;
329:   PetscFunctionBeginUser;
330:   /* get coords for the element */
331:   for (i = 0; i < 4; i++) {
332:     for (d = 0; d < NSD; d++) el_coords[NSD * i + d] = _coords[NSD * e2n[i] + d];
333:   }
334:   PetscFunctionReturn(PETSC_SUCCESS);
335: }

337: static PetscErrorCode AssembleStokes_A(Mat A, DM stokes_da, DM quadrature)
338: {
339:   DM                 cda;
340:   Vec                coords;
341:   const PetscScalar *_coords;
342:   PetscInt           u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
343:   PetscInt           p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
344:   PetscInt           nel, npe, eidx;
345:   const PetscInt    *element_list;
346:   PetscScalar        Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
347:   PetscScalar        Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
348:   PetscScalar        De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
349:   PetscScalar        Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
350:   PetscScalar        el_coords[NODES_PER_EL * NSD];
351:   PetscScalar       *q_eta, *prop_eta;

353:   PetscFunctionBeginUser;
354:   PetscCall(MatZeroEntries(A));
355:   /* setup for coords */
356:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
357:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
358:   PetscCall(VecGetArrayRead(coords, &_coords));

360:   /* setup for coefficients */
361:   PetscCall(DMSwarmGetField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));

363:   PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
364:   for (eidx = 0; eidx < nel; eidx++) {
365:     const PetscInt *element = &element_list[npe * eidx];

367:     /* get coords for the element */
368:     PetscCall(GetElementCoords(_coords, element, el_coords));

370:     /* get coefficients for the element */
371:     prop_eta = &q_eta[GAUSS_POINTS * eidx];

373:     /* initialise element stiffness matrix */
374:     PetscCall(PetscMemzero(Ae, sizeof(Ae)));
375:     PetscCall(PetscMemzero(Ge, sizeof(Ge)));
376:     PetscCall(PetscMemzero(De, sizeof(De)));
377:     PetscCall(PetscMemzero(Ce, sizeof(Ce)));

379:     /* form element stiffness matrix */
380:     BForm_DivT(Ae, el_coords, prop_eta);
381:     BForm_Grad(Ge, el_coords);
382:     BForm_Div(De, el_coords);
383:     BForm_Stabilisation(Ce, el_coords, prop_eta);

385:     /* insert element matrix into global matrix */
386:     PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
387:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
388:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
389:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
390:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
391:   }
392:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
393:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

395:   PetscCall(DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
396:   PetscCall(VecRestoreArrayRead(coords, &_coords));
397:   PetscFunctionReturn(PETSC_SUCCESS);
398: }

400: static PetscErrorCode AssembleStokes_PC(Mat A, DM stokes_da, DM quadrature)
401: {
402:   DM                 cda;
403:   Vec                coords;
404:   const PetscScalar *_coords;
405:   PetscInt           u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
406:   PetscInt           p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
407:   PetscInt           nel, npe, eidx;
408:   const PetscInt    *element_list;
409:   PetscScalar        Ae[NODES_PER_EL * U_DOFS * NODES_PER_EL * U_DOFS];
410:   PetscScalar        Ge[NODES_PER_EL * U_DOFS * NODES_PER_EL * P_DOFS];
411:   PetscScalar        De[NODES_PER_EL * P_DOFS * NODES_PER_EL * U_DOFS];
412:   PetscScalar        Ce[NODES_PER_EL * P_DOFS * NODES_PER_EL * P_DOFS];
413:   PetscScalar        el_coords[NODES_PER_EL * NSD];
414:   PetscScalar       *q_eta, *prop_eta;

416:   PetscFunctionBeginUser;
417:   PetscCall(MatZeroEntries(A));
418:   /* setup for coords */
419:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
420:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
421:   PetscCall(VecGetArrayRead(coords, &_coords));

423:   /* setup for coefficients */
424:   PetscCall(DMSwarmGetField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));

426:   PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
427:   for (eidx = 0; eidx < nel; eidx++) {
428:     const PetscInt *element = &element_list[npe * eidx];

430:     /* get coords for the element */
431:     PetscCall(GetElementCoords(_coords, element, el_coords));

433:     /* get coefficients for the element */
434:     prop_eta = &q_eta[GAUSS_POINTS * eidx];

436:     /* initialise element stiffness matrix */
437:     PetscCall(PetscMemzero(Ae, sizeof(Ae)));
438:     PetscCall(PetscMemzero(Ge, sizeof(Ge)));
439:     PetscCall(PetscMemzero(De, sizeof(De)));
440:     PetscCall(PetscMemzero(Ce, sizeof(Ce)));

442:     /* form element stiffness matrix */
443:     BForm_DivT(Ae, el_coords, prop_eta);
444:     BForm_Grad(Ge, el_coords);
445:     BForm_ScaledMassMatrix(Ce, el_coords, prop_eta);

447:     /* insert element matrix into global matrix */
448:     PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));
449:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES));
450:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES));
451:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES));
452:     PetscCall(MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES));
453:   }
454:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
455:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

457:   PetscCall(DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
458:   PetscCall(VecRestoreArrayRead(coords, &_coords));

460:   PetscFunctionReturn(PETSC_SUCCESS);
461: }

463: static PetscErrorCode AssembleStokes_RHS(Vec F, DM stokes_da, DM quadrature)
464: {
465:   DM                 cda;
466:   Vec                coords;
467:   const PetscScalar *_coords;
468:   PetscInt           u_eqn[NODES_PER_EL * U_DOFS]; /* 2 degrees of freedom */
469:   PetscInt           p_eqn[NODES_PER_EL * P_DOFS]; /* 1 degrees of freedom */
470:   PetscInt           nel, npe, eidx, i;
471:   const PetscInt    *element_list;
472:   PetscScalar        Fe[NODES_PER_EL * U_DOFS];
473:   PetscScalar        He[NODES_PER_EL * P_DOFS];
474:   PetscScalar        el_coords[NODES_PER_EL * NSD];
475:   PetscScalar       *q_rhs, *prop_fy;
476:   Vec                local_F;
477:   PetscScalar       *LA_F;

479:   PetscFunctionBeginUser;
480:   PetscCall(VecZeroEntries(F));
481:   /* setup for coords */
482:   PetscCall(DMGetCoordinateDM(stokes_da, &cda));
483:   PetscCall(DMGetCoordinatesLocal(stokes_da, &coords));
484:   PetscCall(VecGetArrayRead(coords, &_coords));

486:   /* setup for coefficients */
487:   PetscCall(DMSwarmGetField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));

489:   /* get access to the vector */
490:   PetscCall(DMGetLocalVector(stokes_da, &local_F));
491:   PetscCall(VecZeroEntries(local_F));
492:   PetscCall(VecGetArray(local_F, &LA_F));

494:   PetscCall(DMDAGetElements(stokes_da, &nel, &npe, &element_list));
495:   for (eidx = 0; eidx < nel; eidx++) {
496:     const PetscInt *element = &element_list[npe * eidx];

498:     /* get coords for the element */
499:     PetscCall(GetElementCoords(_coords, element, el_coords));

501:     /* get coefficients for the element */
502:     prop_fy = &q_rhs[GAUSS_POINTS * eidx];

504:     /* initialise element stiffness matrix */
505:     PetscCall(PetscMemzero(Fe, sizeof(Fe)));
506:     PetscCall(PetscMemzero(He, sizeof(He)));

508:     /* form element stiffness matrix */
509:     LForm_MomentumRHS(Fe, el_coords, NULL, prop_fy);

511:     /* insert element matrix into global matrix */
512:     PetscCall(DMDAGetElementEqnums_up(element, u_eqn, p_eqn));

514:     for (i = 0; i < NODES_PER_EL * U_DOFS; i++) LA_F[u_eqn[i]] += Fe[i];
515:   }
516:   PetscCall(DMSwarmRestoreField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
517:   PetscCall(VecRestoreArrayRead(coords, &_coords));

519:   PetscCall(VecRestoreArray(local_F, &LA_F));
520:   PetscCall(DMLocalToGlobalBegin(stokes_da, local_F, ADD_VALUES, F));
521:   PetscCall(DMLocalToGlobalEnd(stokes_da, local_F, ADD_VALUES, F));
522:   PetscCall(DMRestoreLocalVector(stokes_da, &local_F));

524:   PetscFunctionReturn(PETSC_SUCCESS);
525: }

527: PetscErrorCode DMSwarmPICInsertPointsCellwise(DM dm, DM dmc, PetscInt e, PetscInt npoints, PetscReal xi[], PetscBool proximity_initialization)
528: {
529:   PetscInt           dim, nel, npe, q, k, d, ncurr;
530:   const PetscInt    *element_list;
531:   Vec                coor;
532:   const PetscScalar *_coor;
533:   PetscReal        **basis, *elcoor, *xp;
534:   PetscReal         *swarm_coor;
535:   PetscInt          *swarm_cellid;

537:   PetscFunctionBeginUser;
538:   PetscCall(DMGetDimension(dm, &dim));
539:   PetscCall(DMDAGetElements(dmc, &nel, &npe, &element_list));

541:   PetscCall(PetscMalloc1(dim * npoints, &xp));
542:   PetscCall(PetscMalloc1(dim * npe, &elcoor));
543:   PetscCall(PetscMalloc1(npoints, &basis));
544:   for (q = 0; q < npoints; q++) {
545:     PetscCall(PetscMalloc1(npe, &basis[q]));

547:     switch (dim) {
548:     case 1:
549:       basis[q][0] = 0.5 * (1.0 - xi[dim * q + 0]);
550:       basis[q][1] = 0.5 * (1.0 + xi[dim * q + 0]);
551:       break;
552:     case 2:
553:       basis[q][0] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
554:       basis[q][1] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]);
555:       basis[q][2] = 0.25 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
556:       basis[q][3] = 0.25 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]);
557:       break;

559:     case 3:
560:       basis[q][0] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
561:       basis[q][1] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
562:       basis[q][2] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
563:       basis[q][3] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 - xi[dim * q + 2]);
564:       basis[q][4] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
565:       basis[q][5] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 - xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
566:       basis[q][6] = 0.125 * (1.0 + xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
567:       basis[q][7] = 0.125 * (1.0 - xi[dim * q + 0]) * (1.0 + xi[dim * q + 1]) * (1.0 + xi[dim * q + 2]);
568:       break;
569:     }
570:   }

572:   PetscCall(DMGetCoordinatesLocal(dmc, &coor));
573:   PetscCall(VecGetArrayRead(coor, &_coor));
574:   /* compute and store the coordinates for the new points */
575:   {
576:     const PetscInt *element = &element_list[npe * e];

578:     for (k = 0; k < npe; k++) {
579:       for (d = 0; d < dim; d++) elcoor[dim * k + d] = PetscRealPart(_coor[dim * element[k] + d]);
580:     }
581:     for (q = 0; q < npoints; q++) {
582:       for (d = 0; d < dim; d++) xp[dim * q + d] = 0.0;
583:       for (k = 0; k < npe; k++) {
584:         for (d = 0; d < dim; d++) xp[dim * q + d] += basis[q][k] * elcoor[dim * k + d];
585:       }
586:     }
587:   }
588:   PetscCall(VecRestoreArrayRead(coor, &_coor));
589:   PetscCall(DMDARestoreElements(dmc, &nel, &npe, &element_list));

591:   PetscCall(DMSwarmGetLocalSize(dm, &ncurr));
592:   PetscCall(DMSwarmAddNPoints(dm, npoints));

594:   if (proximity_initialization) {
595:     PetscInt  *nnlist;
596:     PetscReal *coor_q, *coor_qn;
597:     PetscInt   npoints_e, *plist_e;

599:     PetscCall(DMSwarmSortGetPointsPerCell(dm, e, &npoints_e, &plist_e));

601:     PetscCall(PetscMalloc1(npoints, &nnlist));
602:     /* find nearest neighbour points in this cell */
603:     PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
604:     PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
605:     for (q = 0; q < npoints; q++) {
606:       PetscInt  qn, nearest_neighbour = -1;
607:       PetscReal sep, min_sep          = PETSC_MAX_REAL;

609:       coor_q = &xp[dim * q];
610:       for (qn = 0; qn < npoints_e; qn++) {
611:         coor_qn = &swarm_coor[dim * plist_e[qn]];
612:         sep     = 0.0;
613:         for (d = 0; d < dim; d++) sep += (coor_q[d] - coor_qn[d]) * (coor_q[d] - coor_qn[d]);
614:         if (sep < min_sep) {
615:           nearest_neighbour = plist_e[qn];
616:           min_sep           = sep;
617:         }
618:       }
619:       PetscCheck(nearest_neighbour != -1, PETSC_COMM_SELF, PETSC_ERR_USER, "Cell %" PetscInt_FMT " is empty - cannot initialize using nearest neighbours", e);
620:       nnlist[q] = nearest_neighbour;
621:     }
622:     PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
623:     PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));

625:     /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
626:     for (q = 0; q < npoints; q++) PetscCall(DMSwarmCopyPoint(dm, nnlist[q], ncurr + q));
627:     PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
628:     PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
629:     for (q = 0; q < npoints; q++) {
630:       /* set the coordinates */
631:       for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
632:       /* set the cell index */
633:       swarm_cellid[ncurr + q] = e;
634:     }
635:     PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
636:     PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));

638:     PetscCall(PetscFree(plist_e));
639:     PetscCall(PetscFree(nnlist));
640:   } else {
641:     PetscCall(DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
642:     PetscCall(DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
643:     for (q = 0; q < npoints; q++) {
644:       /* set the coordinates */
645:       for (d = 0; d < dim; d++) swarm_coor[dim * (ncurr + q) + d] = xp[dim * q + d];
646:       /* set the cell index */
647:       swarm_cellid[ncurr + q] = e;
648:     }
649:     PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid));
650:     PetscCall(DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor));
651:   }

653:   PetscCall(PetscFree(xp));
654:   PetscCall(PetscFree(elcoor));
655:   for (q = 0; q < npoints; q++) PetscCall(PetscFree(basis[q]));
656:   PetscCall(PetscFree(basis));
657:   PetscFunctionReturn(PETSC_SUCCESS);
658: }

660: PetscErrorCode MaterialPoint_PopulateCell(DM dm_vp, DM dm_mpoint)
661: {
662:   PetscInt         _npe, _nel, e, nel;
663:   const PetscInt  *element;
664:   DM               dmc;
665:   PetscQuadrature  quadrature;
666:   const PetscReal *xi;
667:   PetscInt         npoints_q, cnt, cnt_g;

669:   PetscFunctionBeginUser;
670:   PetscCall(DMDAGetElements(dm_vp, &_nel, &_npe, &element));
671:   nel = _nel;
672:   PetscCall(DMDARestoreElements(dm_vp, &_nel, &_npe, &element));

674:   PetscCall(PetscDTGaussTensorQuadrature(2, 1, 4, -1.0, 1.0, &quadrature));
675:   PetscCall(PetscQuadratureGetData(quadrature, NULL, NULL, &npoints_q, &xi, NULL));
676:   PetscCall(DMSwarmGetCellDM(dm_mpoint, &dmc));

678:   PetscCall(DMSwarmSortGetAccess(dm_mpoint));

680:   cnt = 0;
681:   for (e = 0; e < nel; e++) {
682:     PetscInt npoints_per_cell;

684:     PetscCall(DMSwarmSortGetNumberOfPointsPerCell(dm_mpoint, e, &npoints_per_cell));

686:     if (npoints_per_cell < 12) {
687:       PetscCall(DMSwarmPICInsertPointsCellwise(dm_mpoint, dm_vp, e, npoints_q, (PetscReal *)xi, PETSC_TRUE));
688:       cnt++;
689:     }
690:   }
691:   PetscCallMPI(MPI_Allreduce(&cnt, &cnt_g, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
692:   if (cnt_g > 0) PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... ....pop cont: adjusted %" PetscInt_FMT " cells\n", cnt_g));

694:   PetscCall(DMSwarmSortRestoreAccess(dm_mpoint));
695:   PetscCall(PetscQuadratureDestroy(&quadrature));
696:   PetscFunctionReturn(PETSC_SUCCESS);
697: }

699: PetscErrorCode MaterialPoint_AdvectRK1(DM dm_vp, Vec vp, PetscReal dt, DM dm_mpoint)
700: {
701:   Vec                vp_l, coor_l;
702:   const PetscScalar *LA_vp;
703:   PetscInt           i, p, e, npoints, nel, npe;
704:   PetscInt          *mpfield_cell;
705:   PetscReal         *mpfield_coor;
706:   const PetscInt    *element_list;
707:   const PetscInt    *element;
708:   PetscScalar        xi_p[NSD], Ni[NODES_PER_EL];
709:   const PetscScalar *LA_coor;
710:   PetscScalar        dx[NSD];

712:   PetscFunctionBeginUser;
713:   PetscCall(DMGetCoordinatesLocal(dm_vp, &coor_l));
714:   PetscCall(VecGetArrayRead(coor_l, &LA_coor));

716:   PetscCall(DMGetLocalVector(dm_vp, &vp_l));
717:   PetscCall(DMGlobalToLocalBegin(dm_vp, vp, INSERT_VALUES, vp_l));
718:   PetscCall(DMGlobalToLocalEnd(dm_vp, vp, INSERT_VALUES, vp_l));
719:   PetscCall(VecGetArrayRead(vp_l, &LA_vp));

721:   PetscCall(DMDAGetElements(dm_vp, &nel, &npe, &element_list));
722:   PetscCall(DMSwarmGetLocalSize(dm_mpoint, &npoints));
723:   PetscCall(DMSwarmGetField(dm_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
724:   PetscCall(DMSwarmGetField(dm_mpoint, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
725:   for (p = 0; p < npoints; p++) {
726:     PetscReal         *coor_p;
727:     PetscScalar        vel_n[NSD * NODES_PER_EL], vel_p[NSD];
728:     const PetscScalar *x0;
729:     const PetscScalar *x2;

731:     e       = mpfield_cell[p];
732:     coor_p  = &mpfield_coor[NSD * p];
733:     element = &element_list[NODES_PER_EL * e];

735:     /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
736:     x0 = &LA_coor[NSD * element[0]];
737:     x2 = &LA_coor[NSD * element[2]];

739:     dx[0] = x2[0] - x0[0];
740:     dx[1] = x2[1] - x0[1];

742:     xi_p[0] = 2.0 * (coor_p[0] - x0[0]) / dx[0] - 1.0;
743:     xi_p[1] = 2.0 * (coor_p[1] - x0[1]) / dx[1] - 1.0;
744:     PetscCheck(PetscRealPart(xi_p[0]) >= -1.0 - PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (xi) too small %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[0]), e);
745:     PetscCheck(PetscRealPart(xi_p[0]) <= 1.0 + PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (xi) too large %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[0]), e);
746:     PetscCheck(PetscRealPart(xi_p[1]) >= -1.0 - PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (eta) too small %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[1]), e);
747:     PetscCheck(PetscRealPart(xi_p[1]) <= 1.0 + PETSC_SMALL, PETSC_COMM_SELF, PETSC_ERR_SUP, "value (eta) too large %1.4e [e=%" PetscInt_FMT "]", (double)PetscRealPart(xi_p[1]), e);

749:     /* evaluate basis functions */
750:     EvaluateBasis_Q1(xi_p, Ni);

752:     /* get cell nodal velocities */
753:     for (i = 0; i < NODES_PER_EL; i++) {
754:       PetscInt nid;

756:       nid                = element[i];
757:       vel_n[NSD * i + 0] = LA_vp[(NSD + 1) * nid + 0];
758:       vel_n[NSD * i + 1] = LA_vp[(NSD + 1) * nid + 1];
759:     }

761:     /* interpolate velocity */
762:     vel_p[0] = vel_p[1] = 0.0;
763:     for (i = 0; i < NODES_PER_EL; i++) {
764:       vel_p[0] += Ni[i] * vel_n[NSD * i + 0];
765:       vel_p[1] += Ni[i] * vel_n[NSD * i + 1];
766:     }

768:     coor_p[0] += dt * PetscRealPart(vel_p[0]);
769:     coor_p[1] += dt * PetscRealPart(vel_p[1]);
770:   }

772:   PetscCall(DMSwarmRestoreField(dm_mpoint, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell));
773:   PetscCall(DMSwarmRestoreField(dm_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor));
774:   PetscCall(DMDARestoreElements(dm_vp, &nel, &npe, &element_list));
775:   PetscCall(VecRestoreArrayRead(vp_l, &LA_vp));
776:   PetscCall(DMRestoreLocalVector(dm_vp, &vp_l));
777:   PetscCall(VecRestoreArrayRead(coor_l, &LA_coor));
778:   PetscFunctionReturn(PETSC_SUCCESS);
779: }

781: PetscErrorCode MaterialPoint_Interpolate(DM dm, Vec eta_v, Vec rho_v, DM dm_quadrature)
782: {
783:   Vec             eta_l, rho_l;
784:   PetscScalar    *_eta_l, *_rho_l;
785:   PetscInt        nqp, npe, nel;
786:   PetscScalar     qp_xi[GAUSS_POINTS][NSD];
787:   PetscScalar     qp_weight[GAUSS_POINTS];
788:   PetscInt        q, k, e;
789:   PetscScalar     Ni[GAUSS_POINTS][NODES_PER_EL];
790:   const PetscInt *element_list;
791:   PetscReal      *q_eta, *q_rhs;

793:   PetscFunctionBeginUser;
794:   /* define quadrature rule */
795:   CreateGaussQuadrature(&nqp, qp_xi, qp_weight);
796:   for (q = 0; q < nqp; q++) EvaluateBasis_Q1(qp_xi[q], Ni[q]);

798:   PetscCall(DMGetLocalVector(dm, &eta_l));
799:   PetscCall(DMGetLocalVector(dm, &rho_l));

801:   PetscCall(DMGlobalToLocalBegin(dm, eta_v, INSERT_VALUES, eta_l));
802:   PetscCall(DMGlobalToLocalEnd(dm, eta_v, INSERT_VALUES, eta_l));
803:   PetscCall(DMGlobalToLocalBegin(dm, rho_v, INSERT_VALUES, rho_l));
804:   PetscCall(DMGlobalToLocalEnd(dm, rho_v, INSERT_VALUES, rho_l));

806:   PetscCall(VecGetArray(eta_l, &_eta_l));
807:   PetscCall(VecGetArray(rho_l, &_rho_l));

809:   PetscCall(DMSwarmGetField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));
810:   PetscCall(DMSwarmGetField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));

812:   PetscCall(DMDAGetElements(dm, &nel, &npe, &element_list));
813:   for (e = 0; e < nel; e++) {
814:     PetscScalar     eta_field_e[NODES_PER_EL];
815:     PetscScalar     rho_field_e[NODES_PER_EL];
816:     const PetscInt *element = &element_list[4 * e];

818:     for (k = 0; k < NODES_PER_EL; k++) {
819:       eta_field_e[k] = _eta_l[element[k]];
820:       rho_field_e[k] = _rho_l[element[k]];
821:     }

823:     for (q = 0; q < nqp; q++) {
824:       PetscScalar eta_q, rho_q;

826:       eta_q = rho_q = 0.0;
827:       for (k = 0; k < NODES_PER_EL; k++) {
828:         eta_q += Ni[q][k] * eta_field_e[k];
829:         rho_q += Ni[q][k] * rho_field_e[k];
830:       }

832:       q_eta[nqp * e + q] = PetscRealPart(eta_q);
833:       q_rhs[nqp * e + q] = PetscRealPart(rho_q);
834:     }
835:   }
836:   PetscCall(DMDARestoreElements(dm, &nel, &npe, &element_list));

838:   PetscCall(DMSwarmRestoreField(dm_quadrature, "rho_q", NULL, NULL, (void **)&q_rhs));
839:   PetscCall(DMSwarmRestoreField(dm_quadrature, "eta_q", NULL, NULL, (void **)&q_eta));

841:   PetscCall(VecRestoreArray(rho_l, &_rho_l));
842:   PetscCall(VecRestoreArray(eta_l, &_eta_l));
843:   PetscCall(DMRestoreLocalVector(dm, &rho_l));
844:   PetscCall(DMRestoreLocalVector(dm, &eta_l));
845:   PetscFunctionReturn(PETSC_SUCCESS);
846: }

848: static PetscErrorCode SolveTimeDepStokes(PetscInt mx, PetscInt my)
849: {
850:   DM              dm_stokes, dm_coeff;
851:   PetscInt        u_dof, p_dof, dof, stencil_width;
852:   Mat             A, B;
853:   PetscInt        nel_local;
854:   Vec             eta_v, rho_v;
855:   Vec             f, X;
856:   KSP             ksp;
857:   PC              pc;
858:   char            filename[PETSC_MAX_PATH_LEN];
859:   DM              dms_quadrature, dms_mpoint;
860:   PetscInt        nel, npe, npoints;
861:   const PetscInt *element_list;
862:   PetscInt        tk, nt, dump_freq;
863:   PetscReal       dt, dt_max = 0.0;
864:   PetscReal       vx[2], vy[2], max_v = 0.0, max_v_step, dh;
865:   const char     *fieldnames[] = {"eta", "rho"};
866:   Vec            *pfields;
867:   PetscInt        ppcell = 1;
868:   PetscReal       time, delta_eta = 1.0;
869:   PetscBool       randomize_coords = PETSC_FALSE;
870:   PetscReal       randomize_fac    = 0.25;
871:   PetscBool       no_view          = PETSC_FALSE;
872:   PetscBool       isbddc;

874:   PetscFunctionBeginUser;
875:   /*
876:     Generate the DMDA for the velocity and pressure spaces.
877:     We use Q1 elements for both fields.
878:     The Q1 FE basis on a regular mesh has a 9-point stencil (DMDA_STENCIL_BOX)
879:     The number of nodes in each direction is mx+1, my+1
880:   */
881:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
882:   p_dof         = P_DOFS; /* p - pressure */
883:   dof           = u_dof + p_dof;
884:   stencil_width = 1;
885:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, mx + 1, my + 1, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, &dm_stokes));
886:   PetscCall(DMDASetElementType(dm_stokes, DMDA_ELEMENT_Q1));
887:   PetscCall(DMSetMatType(dm_stokes, MATAIJ));
888:   PetscCall(DMSetFromOptions(dm_stokes));
889:   PetscCall(DMSetUp(dm_stokes));
890:   PetscCall(DMDASetFieldName(dm_stokes, 0, "ux"));
891:   PetscCall(DMDASetFieldName(dm_stokes, 1, "uy"));
892:   PetscCall(DMDASetFieldName(dm_stokes, 2, "p"));

894:   /* unit box [0,0.9142] x [0,1] */
895:   PetscCall(DMDASetUniformCoordinates(dm_stokes, 0.0, 0.9142, 0.0, 1.0, 0., 0.));
896:   dh = 1.0 / ((PetscReal)(mx));

898:   /* Get local number of elements */
899:   {
900:     PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));

902:     nel_local = nel;

904:     PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
905:   }

907:   /* Create DMDA for representing scalar fields */
908:   PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 1, &dm_coeff));

910:   /* Create the swarm for storing quadrature point values */
911:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_quadrature));
912:   PetscCall(DMSetType(dms_quadrature, DMSWARM));
913:   PetscCall(DMSetDimension(dms_quadrature, 2));
914:   PetscCall(PetscObjectSetName((PetscObject)dms_quadrature, "Quadrature Swarm"));

916:   /* Register fields for viscosity and density on the quadrature points */
917:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "eta_q", 1, PETSC_REAL));
918:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms_quadrature, "rho_q", 1, PETSC_REAL));
919:   PetscCall(DMSwarmFinalizeFieldRegister(dms_quadrature));
920:   PetscCall(DMSwarmSetLocalSizes(dms_quadrature, nel_local * GAUSS_POINTS, 0));

922:   /* Create the material point swarm */
923:   PetscCall(DMCreate(PETSC_COMM_WORLD, &dms_mpoint));
924:   PetscCall(DMSetType(dms_mpoint, DMSWARM));
925:   PetscCall(DMSetDimension(dms_mpoint, 2));
926:   PetscCall(PetscObjectSetName((PetscObject)dms_mpoint, "Material Point Swarm"));

928:   /* Configure the material point swarm to be of type Particle-In-Cell */
929:   PetscCall(DMSwarmSetType(dms_mpoint, DMSWARM_PIC));

931:   /*
932:      Specify the DM to use for point location and projections
933:      within the context of a PIC scheme
934:   */
935:   PetscCall(DMSwarmSetCellDM(dms_mpoint, dm_coeff));

937:   /* Register fields for viscosity and density */
938:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "eta", 1, PETSC_REAL));
939:   PetscCall(DMSwarmRegisterPetscDatatypeField(dms_mpoint, "rho", 1, PETSC_REAL));
940:   PetscCall(DMSwarmFinalizeFieldRegister(dms_mpoint));

942:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL));
943:   PetscCall(DMSwarmSetLocalSizes(dms_mpoint, nel_local * ppcell, 100));

945:   /*
946:     Layout the material points in space using the cell DM.
947:     Particle coordinates are defined by cell wise using different methods.
948:     - DMSWARMPIC_LAYOUT_GAUSS defines particles coordinates at the positions
949:                               corresponding to a Gauss quadrature rule with
950:                               ppcell points in each direction.
951:     - DMSWARMPIC_LAYOUT_REGULAR defines particle coordinates at the centoid of
952:                                 ppcell x ppcell quadralaterals defined within the
953:                                 reference element.
954:     - DMSWARMPIC_LAYOUT_SUBDIVISION defines particles coordinates at the centroid
955:                                     of each quadralateral obtained by sub-dividing
956:                                     the reference element cell ppcell times.
957:   */
958:   PetscCall(DMSwarmInsertPointsUsingCellDM(dms_mpoint, DMSWARMPIC_LAYOUT_SUBDIVISION, ppcell));

960:   /*
961:     Defne a high resolution layer of material points across the material interface
962:   */
963:   {
964:     PetscInt  npoints_dir_x[2];
965:     PetscReal min[2], max[2];

967:     npoints_dir_x[0] = (PetscInt)(0.9142 / (0.05 * dh));
968:     npoints_dir_x[1] = (PetscInt)((0.25 - 0.15) / (0.05 * dh));
969:     min[0]           = 0.0;
970:     max[0]           = 0.9142;
971:     min[1]           = 0.05;
972:     max[1]           = 0.35;
973:     PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
974:   }

976:   /*
977:     Define a high resolution layer of material points near the surface of the domain
978:     to deal with weakly compressible Q1-Q1 elements. These elements "self compact"
979:     when applied to buouyancy driven flow. The error in div(u) is O(h).
980:   */
981:   {
982:     PetscInt  npoints_dir_x[2];
983:     PetscReal min[2], max[2];

985:     npoints_dir_x[0] = (PetscInt)(0.9142 / (0.25 * dh));
986:     npoints_dir_x[1] = (PetscInt)(3.0 * dh / (0.25 * dh));
987:     min[0]           = 0.0;
988:     max[0]           = 0.9142;
989:     min[1]           = 1.0 - 3.0 * dh;
990:     max[1]           = 1.0 - 0.0001;
991:     PetscCall(DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES));
992:   }

994:   PetscCall(DMView(dms_mpoint, PETSC_VIEWER_STDOUT_WORLD));

996:   /* Define initial material properties on each particle in the material point swarm */
997:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-delta_eta", &delta_eta, NULL));
998:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-randomize_coords", &randomize_coords, NULL));
999:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-randomize_fac", &randomize_fac, NULL));
1000:   PetscCheck(randomize_fac <= 1.0, PETSC_COMM_WORLD, PETSC_ERR_USER, "The value of -randomize_fac should be <= 1.0");
1001:   {
1002:     PetscReal  *array_x, *array_e, *array_r;
1003:     PetscInt    p;
1004:     PetscRandom r;
1005:     PetscMPIInt rank;

1007:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

1009:     PetscCall(PetscRandomCreate(PETSC_COMM_SELF, &r));
1010:     PetscCall(PetscRandomSetInterval(r, -randomize_fac * dh, randomize_fac * dh));
1011:     PetscCall(PetscRandomSetSeed(r, (unsigned long)rank));
1012:     PetscCall(PetscRandomSeed(r));

1014:     PetscCall(DMDAGetElements(dm_stokes, &nel, &npe, &element_list));

1016:     /*
1017:        Fetch the registered data from the material point DMSwarm.
1018:        The fields "eta" and "rho" were registered by this example.
1019:        The field identified by the the variable DMSwarmPICField_coor
1020:        was registered by the DMSwarm implementation when the function
1021:          DMSwarmSetType(dms_mpoint,DMSWARM_PIC)
1022:        was called. The returned array defines the coordinates of each
1023:        material point in the point swarm.
1024:     */
1025:     PetscCall(DMSwarmGetField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));
1026:     PetscCall(DMSwarmGetField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
1027:     PetscCall(DMSwarmGetField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));

1029:     PetscCall(DMSwarmGetLocalSize(dms_mpoint, &npoints));
1030:     for (p = 0; p < npoints; p++) {
1031:       PetscReal x_p[2], rr[2];

1033:       if (randomize_coords) {
1034:         PetscCall(PetscRandomGetValueReal(r, &rr[0]));
1035:         PetscCall(PetscRandomGetValueReal(r, &rr[1]));
1036:         array_x[2 * p + 0] += rr[0];
1037:         array_x[2 * p + 1] += rr[1];
1038:       }

1040:       /* Get the coordinates of point, p */
1041:       x_p[0] = array_x[2 * p + 0];
1042:       x_p[1] = array_x[2 * p + 1];

1044:       if (x_p[1] < (0.2 + 0.02 * PetscCosReal(PETSC_PI * x_p[0] / 0.9142))) {
1045:         /* Material properties below the interface */
1046:         array_e[p] = 1.0 * (1.0 / delta_eta);
1047:         array_r[p] = 0.0;
1048:       } else {
1049:         /* Material properties above the interface */
1050:         array_e[p] = 1.0;
1051:         array_r[p] = 1.0;
1052:       }
1053:     }

1055:     /*
1056:        Restore the fetched data fields from the material point DMSwarm.
1057:        Calling the Restore function invalidates the points array_r, array_e, array_x
1058:        by setting them to NULL.
1059:     */
1060:     PetscCall(DMSwarmRestoreField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r));
1061:     PetscCall(DMSwarmRestoreField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e));
1062:     PetscCall(DMSwarmRestoreField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x));

1064:     PetscCall(DMDARestoreElements(dm_stokes, &nel, &npe, &element_list));
1065:     PetscCall(PetscRandomDestroy(&r));
1066:   }

1068:   /*
1069:      If the particle coordinates where randomly shifted, they may have crossed into another
1070:      element, or into another sub-domain. To account for this we call the Migrate function.
1071:   */
1072:   if (randomize_coords) PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));

1074:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-no_view", &no_view, NULL));
1075:   if (!no_view) PetscCall(DMSwarmViewXDMF(dms_mpoint, "ic_coeff_dms.xmf"));

1077:   /* project the swarm properties */
1078:   PetscCall(DMSwarmProjectFields(dms_mpoint, 2, fieldnames, &pfields, PETSC_FALSE));
1079:   eta_v = pfields[0];
1080:   rho_v = pfields[1];
1081:   PetscCall(PetscObjectSetName((PetscObject)eta_v, "eta"));
1082:   PetscCall(PetscObjectSetName((PetscObject)rho_v, "rho"));
1083:   PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));

1085:   /* view projected coefficients eta and rho */
1086:   if (!no_view) {
1087:     PetscViewer viewer;

1089:     PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1090:     PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1091:     PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1092:     PetscCall(PetscViewerFileSetName(viewer, "ic_coeff_dmda.vts"));
1093:     PetscCall(VecView(eta_v, viewer));
1094:     PetscCall(VecView(rho_v, viewer));
1095:     PetscCall(PetscViewerDestroy(&viewer));
1096:   }

1098:   PetscCall(DMCreateMatrix(dm_stokes, &A));
1099:   PetscCall(DMCreateMatrix(dm_stokes, &B));
1100:   PetscCall(DMCreateGlobalVector(dm_stokes, &f));
1101:   PetscCall(DMCreateGlobalVector(dm_stokes, &X));

1103:   PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
1104:   PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
1105:   PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));

1107:   PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
1108:   PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));

1110:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
1111:   PetscCall(KSPSetOptionsPrefix(ksp, "stokes_"));
1112:   PetscCall(KSPSetDM(ksp, dm_stokes));
1113:   PetscCall(KSPSetDMActive(ksp, PETSC_FALSE));
1114:   PetscCall(KSPSetOperators(ksp, A, B));
1115:   PetscCall(KSPSetFromOptions(ksp));
1116:   PetscCall(KSPGetPC(ksp, &pc));
1117:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &isbddc));
1118:   if (isbddc) PetscCall(KSPSetOperators(ksp, A, A));

1120:   /* Define u-v-p indices for fieldsplit */
1121:   {
1122:     PC             pc;
1123:     const PetscInt ufields[] = {0, 1}, pfields[1] = {2};

1125:     PetscCall(KSPGetPC(ksp, &pc));
1126:     PetscCall(PCFieldSplitSetBlockSize(pc, 3));
1127:     PetscCall(PCFieldSplitSetFields(pc, "u", 2, ufields, ufields));
1128:     PetscCall(PCFieldSplitSetFields(pc, "p", 1, pfields, pfields));
1129:   }

1131:   /* If using a fieldsplit preconditioner, attach a DMDA to the velocity split so that geometric multigrid can be used */
1132:   {
1133:     PC        pc, pc_u;
1134:     KSP      *sub_ksp, ksp_u;
1135:     PetscInt  nsplits;
1136:     DM        dm_u;
1137:     PetscBool is_pcfs;

1139:     PetscCall(KSPGetPC(ksp, &pc));

1141:     is_pcfs = PETSC_FALSE;
1142:     PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs));

1144:     if (is_pcfs) {
1145:       PetscCall(KSPSetUp(ksp));
1146:       PetscCall(KSPGetPC(ksp, &pc));
1147:       PetscCall(PCFieldSplitGetSubKSP(pc, &nsplits, &sub_ksp));
1148:       ksp_u = sub_ksp[0];
1149:       PetscCall(PetscFree(sub_ksp));

1151:       if (nsplits == 2) {
1152:         PetscCall(DMDACreateCompatibleDMDA(dm_stokes, 2, &dm_u));

1154:         PetscCall(KSPSetDM(ksp_u, dm_u));
1155:         PetscCall(KSPSetDMActive(ksp_u, PETSC_FALSE));
1156:         PetscCall(DMDestroy(&dm_u));

1158:         /* enforce galerkin coarse grids be used */
1159:         PetscCall(KSPGetPC(ksp_u, &pc_u));
1160:         PetscCall(PCMGSetGalerkin(pc_u, PC_MG_GALERKIN_PMAT));
1161:       }
1162:     }
1163:   }

1165:   dump_freq = 10;
1166:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dump_freq", &dump_freq, NULL));
1167:   nt = 10;
1168:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nt", &nt, NULL));
1169:   time = 0.0;
1170:   for (tk = 1; tk <= nt; tk++) {
1171:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... assemble\n"));
1172:     PetscCall(AssembleStokes_A(A, dm_stokes, dms_quadrature));
1173:     PetscCall(AssembleStokes_PC(B, dm_stokes, dms_quadrature));
1174:     PetscCall(AssembleStokes_RHS(f, dm_stokes, dms_quadrature));

1176:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... bc imposition\n"));
1177:     PetscCall(DMDAApplyBoundaryConditions(dm_stokes, A, f));
1178:     PetscCall(DMDAApplyBoundaryConditions(dm_stokes, B, NULL));

1180:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... solve\n"));
1181:     PetscCall(KSPSetOperators(ksp, A, isbddc ? A : B));
1182:     PetscCall(KSPSolve(ksp, f, X));

1184:     PetscCall(VecStrideMax(X, 0, NULL, &vx[1]));
1185:     PetscCall(VecStrideMax(X, 1, NULL, &vy[1]));
1186:     PetscCall(VecStrideMin(X, 0, NULL, &vx[0]));
1187:     PetscCall(VecStrideMin(X, 1, NULL, &vy[0]));

1189:     max_v_step = PetscMax(vx[0], vx[1]);
1190:     max_v_step = PetscMax(max_v_step, vy[0]);
1191:     max_v_step = PetscMax(max_v_step, vy[1]);
1192:     max_v      = PetscMax(max_v, max_v_step);

1194:     dt_max = 2.0;
1195:     dt     = 0.5 * (dh / max_v_step);
1196:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... max v %1.4e , dt %1.4e : [total] max v %1.4e , dt_max %1.4e\n", (double)max_v_step, (double)dt, (double)max_v, (double)dt_max));
1197:     dt = PetscMin(dt_max, dt);

1199:     /* advect */
1200:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... advect\n"));
1201:     PetscCall(MaterialPoint_AdvectRK1(dm_stokes, X, dt, dms_mpoint));

1203:     /* migrate */
1204:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... migrate\n"));
1205:     PetscCall(DMSwarmMigrate(dms_mpoint, PETSC_TRUE));

1207:     /* update cell population */
1208:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... populate cells\n"));
1209:     PetscCall(MaterialPoint_PopulateCell(dm_stokes, dms_mpoint));

1211:     /* update coefficients on quadrature points */
1212:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... project\n"));
1213:     PetscCall(DMSwarmProjectFields(dms_mpoint, 2, fieldnames, &pfields, PETSC_TRUE));
1214:     eta_v = pfields[0];
1215:     rho_v = pfields[1];
1216:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... interp\n"));
1217:     PetscCall(MaterialPoint_Interpolate(dm_coeff, eta_v, rho_v, dms_quadrature));

1219:     if (tk % dump_freq == 0) {
1220:       PetscViewer viewer;

1222:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, ".... write XDMF, VTS\n"));
1223:       PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_coeff_dms.xmf", tk));
1224:       PetscCall(DMSwarmViewXDMF(dms_mpoint, filename));

1226:       PetscCall(PetscSNPrintf(filename, PETSC_MAX_PATH_LEN - 1, "step%.4" PetscInt_FMT "_vp_dm.vts", tk));
1227:       PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
1228:       PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
1229:       PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
1230:       PetscCall(PetscViewerFileSetName(viewer, filename));
1231:       PetscCall(VecView(X, viewer));
1232:       PetscCall(PetscViewerDestroy(&viewer));
1233:     }
1234:     time += dt;
1235:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step %" PetscInt_FMT " : time %1.2e \n", tk, (double)time));
1236:   }

1238:   PetscCall(KSPDestroy(&ksp));
1239:   PetscCall(VecDestroy(&X));
1240:   PetscCall(VecDestroy(&f));
1241:   PetscCall(MatDestroy(&A));
1242:   PetscCall(MatDestroy(&B));
1243:   PetscCall(VecDestroy(&eta_v));
1244:   PetscCall(VecDestroy(&rho_v));
1245:   PetscCall(PetscFree(pfields));

1247:   PetscCall(DMDestroy(&dms_mpoint));
1248:   PetscCall(DMDestroy(&dms_quadrature));
1249:   PetscCall(DMDestroy(&dm_coeff));
1250:   PetscCall(DMDestroy(&dm_stokes));
1251:   PetscFunctionReturn(PETSC_SUCCESS);
1252: }

1254: /*
1255:  <sequential run>
1256:  ./ex70 -stokes_ksp_type fgmres -stokes_pc_type fieldsplit -stokes_pc_fieldsplit_block_size 3 -stokes_pc_fieldsplit_type SYMMETRIC_MULTIPLICATIVE -stokes_pc_fieldsplit_0_fields 0,1 -stokes_pc_fieldsplit_1_fields 2 -stokes_fieldsplit_0_ksp_type preonly -stokes_fieldsplit_0_pc_type lu -stokes_fieldsplit_1_ksp_type preonly -stokes_fieldsplit_1_pc_type lu  -mx 80 -my 80  -stokes_ksp_converged_reason  -dump_freq 25  -stokes_ksp_rtol 1.0e-8 -build_twosided allreduce  -ppcell 2 -nt 4000 -delta_eta 1.0 -randomize_coords
1257: */
1258: int main(int argc, char **args)
1259: {
1260:   PetscInt  mx, my;
1261:   PetscBool set = PETSC_FALSE;

1263:   PetscFunctionBeginUser;
1264:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
1265:   mx = my = 10;
1266:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL));
1267:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL));
1268:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mxy", &mx, &set));
1269:   if (set) my = mx;
1270:   PetscCall(SolveTimeDepStokes(mx, my));
1271:   PetscCall(PetscFinalize());
1272:   return 0;
1273: }

1275: /* -------------------------- helpers for boundary conditions -------------------------------- */
1276: static PetscErrorCode BCApplyZero_EAST(DM da, PetscInt d_idx, Mat A, Vec b)
1277: {
1278:   DM                     cda;
1279:   Vec                    coords;
1280:   PetscInt               si, sj, nx, ny, i, j;
1281:   PetscInt               M, N;
1282:   DMDACoor2d           **_coords;
1283:   const PetscInt        *g_idx;
1284:   PetscInt              *bc_global_ids;
1285:   PetscScalar           *bc_vals;
1286:   PetscInt               nbcs;
1287:   PetscInt               n_dofs;
1288:   ISLocalToGlobalMapping ltogm;

1290:   PetscFunctionBeginUser;
1291:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
1292:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

1294:   PetscCall(DMGetCoordinateDM(da, &cda));
1295:   PetscCall(DMGetCoordinatesLocal(da, &coords));
1296:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1297:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1298:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

1300:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1301:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

1303:   /* init the entries to -1 so VecSetValues will ignore them */
1304:   for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;

1306:   i = nx - 1;
1307:   for (j = 0; j < ny; j++) {
1308:     PetscInt local_id;

1310:     local_id = i + j * nx;

1312:     bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];

1314:     bc_vals[j] = 0.0;
1315:   }
1316:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1317:   nbcs = 0;
1318:   if ((si + nx) == (M)) nbcs = ny;

1320:   if (b) {
1321:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1322:     PetscCall(VecAssemblyBegin(b));
1323:     PetscCall(VecAssemblyEnd(b));
1324:   }
1325:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

1327:   PetscCall(PetscFree(bc_vals));
1328:   PetscCall(PetscFree(bc_global_ids));

1330:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1331:   PetscFunctionReturn(PETSC_SUCCESS);
1332: }

1334: static PetscErrorCode BCApplyZero_WEST(DM da, PetscInt d_idx, Mat A, Vec b)
1335: {
1336:   DM                     cda;
1337:   Vec                    coords;
1338:   PetscInt               si, sj, nx, ny, i, j;
1339:   PetscInt               M, N;
1340:   DMDACoor2d           **_coords;
1341:   const PetscInt        *g_idx;
1342:   PetscInt              *bc_global_ids;
1343:   PetscScalar           *bc_vals;
1344:   PetscInt               nbcs;
1345:   PetscInt               n_dofs;
1346:   ISLocalToGlobalMapping ltogm;

1348:   PetscFunctionBeginUser;
1349:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
1350:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

1352:   PetscCall(DMGetCoordinateDM(da, &cda));
1353:   PetscCall(DMGetCoordinatesLocal(da, &coords));
1354:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1355:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1356:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

1358:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_global_ids));
1359:   PetscCall(PetscMalloc1(ny * n_dofs, &bc_vals));

1361:   /* init the entries to -1 so VecSetValues will ignore them */
1362:   for (i = 0; i < ny * n_dofs; i++) bc_global_ids[i] = -1;

1364:   i = 0;
1365:   for (j = 0; j < ny; j++) {
1366:     PetscInt local_id;

1368:     local_id = i + j * nx;

1370:     bc_global_ids[j] = g_idx[n_dofs * local_id + d_idx];

1372:     bc_vals[j] = 0.0;
1373:   }
1374:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1375:   nbcs = 0;
1376:   if (si == 0) nbcs = ny;

1378:   if (b) {
1379:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1380:     PetscCall(VecAssemblyBegin(b));
1381:     PetscCall(VecAssemblyEnd(b));
1382:   }

1384:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

1386:   PetscCall(PetscFree(bc_vals));
1387:   PetscCall(PetscFree(bc_global_ids));

1389:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1390:   PetscFunctionReturn(PETSC_SUCCESS);
1391: }

1393: static PetscErrorCode BCApplyZero_NORTH(DM da, PetscInt d_idx, Mat A, Vec b)
1394: {
1395:   DM                     cda;
1396:   Vec                    coords;
1397:   PetscInt               si, sj, nx, ny, i, j;
1398:   PetscInt               M, N;
1399:   DMDACoor2d           **_coords;
1400:   const PetscInt        *g_idx;
1401:   PetscInt              *bc_global_ids;
1402:   PetscScalar           *bc_vals;
1403:   PetscInt               nbcs;
1404:   PetscInt               n_dofs;
1405:   ISLocalToGlobalMapping ltogm;

1407:   PetscFunctionBeginUser;
1408:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
1409:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

1411:   PetscCall(DMGetCoordinateDM(da, &cda));
1412:   PetscCall(DMGetCoordinatesLocal(da, &coords));
1413:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1414:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1415:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

1417:   PetscCall(PetscMalloc1(nx, &bc_global_ids));
1418:   PetscCall(PetscMalloc1(nx, &bc_vals));

1420:   /* init the entries to -1 so VecSetValues will ignore them */
1421:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

1423:   j = ny - 1;
1424:   for (i = 0; i < nx; i++) {
1425:     PetscInt local_id;

1427:     local_id = i + j * nx;

1429:     bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];

1431:     bc_vals[i] = 0.0;
1432:   }
1433:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1434:   nbcs = 0;
1435:   if ((sj + ny) == (N)) nbcs = nx;

1437:   if (b) {
1438:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1439:     PetscCall(VecAssemblyBegin(b));
1440:     PetscCall(VecAssemblyEnd(b));
1441:   }
1442:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, NULL, NULL));

1444:   PetscCall(PetscFree(bc_vals));
1445:   PetscCall(PetscFree(bc_global_ids));

1447:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1448:   PetscFunctionReturn(PETSC_SUCCESS);
1449: }

1451: static PetscErrorCode BCApplyZero_SOUTH(DM da, PetscInt d_idx, Mat A, Vec b)
1452: {
1453:   DM                     cda;
1454:   Vec                    coords;
1455:   PetscInt               si, sj, nx, ny, i, j;
1456:   PetscInt               M, N;
1457:   DMDACoor2d           **_coords;
1458:   const PetscInt        *g_idx;
1459:   PetscInt              *bc_global_ids;
1460:   PetscScalar           *bc_vals;
1461:   PetscInt               nbcs;
1462:   PetscInt               n_dofs;
1463:   ISLocalToGlobalMapping ltogm;

1465:   PetscFunctionBeginUser;
1466:   PetscCall(DMGetLocalToGlobalMapping(da, &ltogm));
1467:   PetscCall(ISLocalToGlobalMappingGetIndices(ltogm, &g_idx));

1469:   PetscCall(DMGetCoordinateDM(da, &cda));
1470:   PetscCall(DMGetCoordinatesLocal(da, &coords));
1471:   PetscCall(DMDAVecGetArray(cda, coords, &_coords));
1472:   PetscCall(DMDAGetGhostCorners(cda, &si, &sj, 0, &nx, &ny, 0));
1473:   PetscCall(DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, &n_dofs, 0, 0, 0, 0, 0));

1475:   PetscCall(PetscMalloc1(nx, &bc_global_ids));
1476:   PetscCall(PetscMalloc1(nx, &bc_vals));

1478:   /* init the entries to -1 so VecSetValues will ignore them */
1479:   for (i = 0; i < nx; i++) bc_global_ids[i] = -1;

1481:   j = 0;
1482:   for (i = 0; i < nx; i++) {
1483:     PetscInt local_id;

1485:     local_id = i + j * nx;

1487:     bc_global_ids[i] = g_idx[n_dofs * local_id + d_idx];

1489:     bc_vals[i] = 0.0;
1490:   }
1491:   PetscCall(ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx));
1492:   nbcs = 0;
1493:   if (sj == 0) nbcs = nx;

1495:   if (b) {
1496:     PetscCall(VecSetValues(b, nbcs, bc_global_ids, bc_vals, INSERT_VALUES));
1497:     PetscCall(VecAssemblyBegin(b));
1498:     PetscCall(VecAssemblyEnd(b));
1499:   }
1500:   if (A) PetscCall(MatZeroRowsColumns(A, nbcs, bc_global_ids, 1.0, 0, 0));

1502:   PetscCall(PetscFree(bc_vals));
1503:   PetscCall(PetscFree(bc_global_ids));

1505:   PetscCall(DMDAVecRestoreArray(cda, coords, &_coords));
1506:   PetscFunctionReturn(PETSC_SUCCESS);
1507: }

1509: /*
1510:  Impose free slip boundary conditions on the left/right faces: u_i n_i = 0, tau_{ij} t_j = 0
1511:  Impose no slip boundray conditions on the top/bottom faces:   u_i n_i = 0, u_i t_i = 0
1512: */
1513: static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes, Mat A, Vec f)
1514: {
1515:   PetscFunctionBeginUser;
1516:   PetscCall(BCApplyZero_NORTH(dm_stokes, 0, A, f));
1517:   PetscCall(BCApplyZero_NORTH(dm_stokes, 1, A, f));
1518:   PetscCall(BCApplyZero_EAST(dm_stokes, 0, A, f));
1519:   PetscCall(BCApplyZero_SOUTH(dm_stokes, 0, A, f));
1520:   PetscCall(BCApplyZero_SOUTH(dm_stokes, 1, A, f));
1521:   PetscCall(BCApplyZero_WEST(dm_stokes, 0, A, f));
1522:   PetscFunctionReturn(PETSC_SUCCESS);
1523: }

1525: /*TEST

1527:    test:
1528:      suffix: 1
1529:      args: -no_view
1530:      requires: !complex double
1531:      filter: grep -v atomic
1532:      filter_output: grep -v atomic
1533:    test:
1534:      suffix: 1_matis
1535:      requires: !complex double
1536:      args: -no_view -dm_mat_type is
1537:      filter: grep -v atomic
1538:      filter_output: grep -v atomic
1539:    testset:
1540:      nsize: 4
1541:      requires: !complex double
1542:      args: -no_view -dm_mat_type is -stokes_ksp_type fetidp -mx 80 -my 80 -stokes_ksp_converged_reason -stokes_ksp_rtol 1.0e-8 -ppcell 2 -nt 4 -randomize_coords -stokes_ksp_error_if_not_converged
1543:      filter: grep -v atomic
1544:      filter_output: grep -v atomic
1545:      test:
1546:        suffix: fetidp
1547:        args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
1548:      test:
1549:        suffix: fetidp_lumped
1550:        args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd -stokes_fetidp_pc_lumped -stokes_fetidp_bddc_pc_bddc_dirichlet_pc_type none -stokes_fetidp_bddc_pc_bddc_switch_static
1551:      test:
1552:        suffix: fetidp_saddlepoint
1553:        args: -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_norm_type natural -stokes_fetidp_pc_fieldsplit_schur_fact_type diag -stokes_fetidp_fieldsplit_p_pc_type bjacobi -stokes_fetidp_fieldsplit_lag_ksp_type preonly -stokes_fetidp_fieldsplit_p_ksp_type preonly -stokes_ksp_fetidp_pressure_field 2 -stokes_fetidp_pc_fieldsplit_schur_scale -1
1554:      test:
1555:        suffix: fetidp_saddlepoint_lumped
1556:        args: -stokes_ksp_fetidp_saddlepoint -stokes_fetidp_ksp_type cg -stokes_ksp_norm_type natural -stokes_fetidp_pc_fieldsplit_schur_fact_type diag -stokes_fetidp_fieldsplit_p_pc_type bjacobi -stokes_fetidp_fieldsplit_lag_ksp_type preonly -stokes_fetidp_fieldsplit_p_ksp_type preonly -stokes_ksp_fetidp_pressure_field 2 -stokes_fetidp_pc_fieldsplit_schur_scale -1 -stokes_fetidp_bddc_pc_bddc_dirichlet_pc_type none -stokes_fetidp_bddc_pc_bddc_switch_static -stokes_fetidp_pc_lumped
1557: TEST*/