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, &ltogm));
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, &ltogm));
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, &ltogm));
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, &ltogm));
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*/