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;
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:   return 0;
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:   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;
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:   return 0;
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;

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

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

363:   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:     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:     PetscMemzero(Ae, sizeof(Ae));
375:     PetscMemzero(Ge, sizeof(Ge));
376:     PetscMemzero(De, sizeof(De));
377:     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:     DMDAGetElementEqnums_up(element, u_eqn, p_eqn);
387:     MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES);
388:     MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES);
389:     MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES);
390:     MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES);
391:   }
392:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
393:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

395:   DMSwarmRestoreField(quadrature, "eta_q", NULL, NULL, (void **)&q_eta);
396:   VecRestoreArrayRead(coords, &_coords);
397:   return 0;
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;

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

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

426:   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:     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:     PetscMemzero(Ae, sizeof(Ae));
438:     PetscMemzero(Ge, sizeof(Ge));
439:     PetscMemzero(De, sizeof(De));
440:     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:     DMDAGetElementEqnums_up(element, u_eqn, p_eqn);
449:     MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * U_DOFS, u_eqn, Ae, ADD_VALUES);
450:     MatSetValuesLocal(A, NODES_PER_EL * U_DOFS, u_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ge, ADD_VALUES);
451:     MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * U_DOFS, u_eqn, De, ADD_VALUES);
452:     MatSetValuesLocal(A, NODES_PER_EL * P_DOFS, p_eqn, NODES_PER_EL * P_DOFS, p_eqn, Ce, ADD_VALUES);
453:   }
454:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
455:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

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

460:   return 0;
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;

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

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

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

494:   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:     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:     PetscMemzero(Fe, sizeof(Fe));
506:     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:     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:   DMSwarmRestoreField(quadrature, "rho_q", NULL, NULL, (void **)&q_rhs);
517:   VecRestoreArrayRead(coords, &_coords);

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

524:   return 0;
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;

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

541:   PetscMalloc1(dim * npoints, &xp);
542:   PetscMalloc1(dim * npe, &elcoor);
543:   PetscMalloc1(npoints, &basis);
544:   for (q = 0; q < npoints; q++) {
545:     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:   DMGetCoordinatesLocal(dmc, &coor);
573:   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:   VecRestoreArrayRead(coor, &_coor);
589:   DMDARestoreElements(dmc, &nel, &npe, &element_list);

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

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

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

601:     PetscMalloc1(npoints, &nnlist);
602:     /* find nearest neighour points in this cell */
603:     DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
604:     DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
605:     for (q = 0; q < npoints; q++) {
606:       PetscInt  qn, nearest_neighour = -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_neighour = plist_e[qn];
616:           min_sep          = sep;
617:         }
618:       }
620:       nnlist[q] = nearest_neighour;
621:     }
622:     DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
623:     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++) DMSwarmCopyPoint(dm, nnlist[q], ncurr + q);
627:     DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
628:     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:     DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
636:     DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);

638:     PetscFree(plist_e);
639:     PetscFree(nnlist);
640:   } else {
641:     DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
642:     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:     DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&swarm_cellid);
650:     DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&swarm_coor);
651:   }

653:   PetscFree(xp);
654:   PetscFree(elcoor);
655:   for (q = 0; q < npoints; q++) PetscFree(basis[q]);
656:   PetscFree(basis);
657:   return 0;
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;

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

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

678:   DMSwarmSortGetAccess(dm_mpoint);

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

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

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

694:   DMSwarmSortRestoreAccess(dm_mpoint);
695:   PetscQuadratureDestroy(&quadrature);
696:   return 0;
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];

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

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

721:   DMDAGetElements(dm_vp, &nel, &npe, &element_list);
722:   DMSwarmGetLocalSize(dm_mpoint, &npoints);
723:   DMSwarmGetField(dm_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor);
724:   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;

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:   DMSwarmRestoreField(dm_mpoint, DMSwarmPICField_cellid, NULL, NULL, (void **)&mpfield_cell);
773:   DMSwarmRestoreField(dm_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&mpfield_coor);
774:   DMDARestoreElements(dm_vp, &nel, &npe, &element_list);
775:   VecRestoreArrayRead(vp_l, &LA_vp);
776:   DMRestoreLocalVector(dm_vp, &vp_l);
777:   VecRestoreArrayRead(coor_l, &LA_coor);
778:   return 0;
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;

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:   DMGetLocalVector(dm, &eta_l);
799:   DMGetLocalVector(dm, &rho_l);

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

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

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

812:   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:   DMDARestoreElements(dm, &nel, &npe, &element_list);

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

841:   VecRestoreArray(rho_l, &_rho_l);
842:   VecRestoreArray(eta_l, &_eta_l);
843:   DMRestoreLocalVector(dm, &rho_l);
844:   DMRestoreLocalVector(dm, &eta_l);
845:   return 0;
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;

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:   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:   DMDASetElementType(dm_stokes, DMDA_ELEMENT_Q1);
887:   DMSetMatType(dm_stokes, MATAIJ);
888:   DMSetFromOptions(dm_stokes);
889:   DMSetUp(dm_stokes);
890:   DMDASetFieldName(dm_stokes, 0, "ux");
891:   DMDASetFieldName(dm_stokes, 1, "uy");
892:   DMDASetFieldName(dm_stokes, 2, "p");

894:   /* unit box [0,0.9142] x [0,1] */
895:   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:     DMDAGetElements(dm_stokes, &nel, &npe, &element_list);

902:     nel_local = nel;

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

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

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

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

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

928:   /* Configure the material point swarm to be of type Particle-In-Cell */
929:   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:   DMSwarmSetCellDM(dms_mpoint, dm_coeff);

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

942:   PetscOptionsGetInt(NULL, NULL, "-ppcell", &ppcell, NULL);
943:   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:   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:     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:     DMSwarmSetPointsUniformCoordinates(dms_mpoint, min, max, npoints_dir_x, ADD_VALUES);
992:   }

994:   DMView(dms_mpoint, PETSC_VIEWER_STDOUT_WORLD);

996:   /* Define initial material properties on each particle in the material point swarm */
997:   PetscOptionsGetReal(NULL, NULL, "-delta_eta", &delta_eta, NULL);
998:   PetscOptionsGetBool(NULL, NULL, "-randomize_coords", &randomize_coords, NULL);
999:   PetscOptionsGetReal(NULL, NULL, "-randomize_fac", &randomize_fac, NULL);
1001:   {
1002:     PetscReal  *array_x, *array_e, *array_r;
1003:     PetscInt    p;
1004:     PetscRandom r;
1005:     PetscMPIInt rank;

1007:     MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

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

1014:     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:     DMSwarmGetField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x);
1026:     DMSwarmGetField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e);
1027:     DMSwarmGetField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r);

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

1033:       if (randomize_coords) {
1034:         PetscRandomGetValueReal(r, &rr[0]);
1035:         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:     DMSwarmRestoreField(dms_mpoint, "rho", NULL, NULL, (void **)&array_r);
1061:     DMSwarmRestoreField(dms_mpoint, "eta", NULL, NULL, (void **)&array_e);
1062:     DMSwarmRestoreField(dms_mpoint, DMSwarmPICField_coor, NULL, NULL, (void **)&array_x);

1064:     DMDARestoreElements(dm_stokes, &nel, &npe, &element_list);
1065:     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) DMSwarmMigrate(dms_mpoint, PETSC_TRUE);

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

1077:   /* project the swarm properties */
1078:   DMSwarmProjectFields(dms_mpoint, 2, fieldnames, &pfields, PETSC_FALSE);
1079:   eta_v = pfields[0];
1080:   rho_v = pfields[1];
1081:   PetscObjectSetName((PetscObject)eta_v, "eta");
1082:   PetscObjectSetName((PetscObject)rho_v, "rho");
1083:   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:     PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
1090:     PetscViewerSetType(viewer, PETSCVIEWERVTK);
1091:     PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
1092:     PetscViewerFileSetName(viewer, "ic_coeff_dmda.vts");
1093:     VecView(eta_v, viewer);
1094:     VecView(rho_v, viewer);
1095:     PetscViewerDestroy(&viewer);
1096:   }

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

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

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

1110:   KSPCreate(PETSC_COMM_WORLD, &ksp);
1111:   KSPSetOptionsPrefix(ksp, "stokes_");
1112:   KSPSetDM(ksp, dm_stokes);
1113:   KSPSetDMActive(ksp, PETSC_FALSE);
1114:   KSPSetOperators(ksp, A, B);
1115:   KSPSetFromOptions(ksp);
1116:   KSPGetPC(ksp, &pc);
1117:   PetscObjectTypeCompare((PetscObject)pc, PCBDDC, &isbddc);
1118:   if (isbddc) 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:     KSPGetPC(ksp, &pc);
1126:     PCFieldSplitSetBlockSize(pc, 3);
1127:     PCFieldSplitSetFields(pc, "u", 2, ufields, ufields);
1128:     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:     KSPGetPC(ksp, &pc);

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

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

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

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

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

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

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

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

1184:     VecStrideMax(X, 0, NULL, &vx[1]);
1185:     VecStrideMax(X, 1, NULL, &vy[1]);
1186:     VecStrideMin(X, 0, NULL, &vx[0]);
1187:     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:     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:     PetscPrintf(PETSC_COMM_WORLD, ".... advect\n");
1201:     MaterialPoint_AdvectRK1(dm_stokes, X, dt, dms_mpoint);

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

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

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

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

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

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

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

1247:   DMDestroy(&dms_mpoint);
1248:   DMDestroy(&dms_quadrature);
1249:   DMDestroy(&dm_coeff);
1250:   DMDestroy(&dm_stokes);
1251:   return 0;
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;

1264:   PetscInitialize(&argc, &args, (char *)0, help);
1265:   mx = my = 10;
1266:   PetscOptionsGetInt(NULL, NULL, "-mx", &mx, NULL);
1267:   PetscOptionsGetInt(NULL, NULL, "-my", &my, NULL);
1268:   PetscOptionsGetInt(NULL, NULL, "-mxy", &mx, &set);
1269:   if (set) my = mx;
1270:   SolveTimeDepStokes(mx, my);
1271:   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;

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

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

1300:   PetscMalloc1(ny * n_dofs, &bc_global_ids);
1301:   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:   ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1317:   nbcs = 0;
1318:   if ((si + nx) == (M)) nbcs = ny;

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

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

1330:   DMDAVecRestoreArray(cda, coords, &_coords);
1331:   return 0;
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;

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

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

1358:   PetscMalloc1(ny * n_dofs, &bc_global_ids);
1359:   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:   ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1375:   nbcs = 0;
1376:   if (si == 0) nbcs = ny;

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

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

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

1389:   DMDAVecRestoreArray(cda, coords, &_coords);
1390:   return 0;
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;

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

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

1417:   PetscMalloc1(nx, &bc_global_ids);
1418:   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:   ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1434:   nbcs = 0;
1435:   if ((sj + ny) == (N)) nbcs = nx;

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

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

1447:   DMDAVecRestoreArray(cda, coords, &_coords);
1448:   return 0;
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;

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

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

1475:   PetscMalloc1(nx, &bc_global_ids);
1476:   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:   ISLocalToGlobalMappingRestoreIndices(ltogm, &g_idx);
1492:   nbcs = 0;
1493:   if (sj == 0) nbcs = nx;

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

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

1505:   DMDAVecRestoreArray(cda, coords, &_coords);
1506:   return 0;
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: {
1516:   BCApplyZero_NORTH(dm_stokes, 0, A, f);
1517:   BCApplyZero_NORTH(dm_stokes, 1, A, f);
1518:   BCApplyZero_EAST(dm_stokes, 0, A, f);
1519:   BCApplyZero_SOUTH(dm_stokes, 0, A, f);
1520:   BCApplyZero_SOUTH(dm_stokes, 1, A, f);
1521:   BCApplyZero_WEST(dm_stokes, 0, A, f);
1522:   return 0;
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*/