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; gp_xi[0][1] = -0.57735026919;
106:   gp_xi[1][0]  = -0.57735026919; gp_xi[1][1] =  0.57735026919;
107:   gp_xi[2][0]  =  0.57735026919; gp_xi[2][1] =  0.57735026919;
108:   gp_xi[3][0]  =  0.57735026919; gp_xi[3][1] = -0.57735026919;
109:   gp_weight[0] = 1.0;
110:   gp_weight[1] = 1.0;
111:   gp_weight[2] = 1.0;
112:   gp_weight[3] = 1.0;
113: }

115: static PetscErrorCode DMDAGetElementEqnums_up(const PetscInt element[],PetscInt s_u[],PetscInt s_p[])
116: {
117:   PetscInt i;
119:   for (i=0; i<NODES_PER_EL; i++) {
120:     /* velocity */
121:     s_u[NSD*i+0] = 3*element[i];
122:     s_u[NSD*i+1] = 3*element[i]+1;
123:     /* pressure */
124:     s_p[i] = 3*element[i]+2;
125:   }
126:   return(0);
127: }

129: 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)
130: {
131:   PetscInt ij,r,c,nc;

133:   nc = u_NPE*u_dof;
134:   r = w_dof*wi+wd;
135:   c = u_dof*ui+ud;
136:   ij = r*nc+c;
137:   return(ij);
138: }

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

148:   /* define quadrature rule */
149:   CreateGaussQuadrature(&ngp,gp_xi,gp_weight);

151:   /* evaluate bilinear form */
152:   for (p = 0; p < ngp; p++) {
153:     EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
154:     EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);

156:     for (i = 0; i < NODES_PER_EL; i++) {
157:       PetscScalar d_dx_i = GNx_p[0][i];
158:       PetscScalar d_dy_i = GNx_p[1][i];

160:       B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
161:       B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
162:       B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
163:     }

165:     tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
166:     tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
167:     tildeD[2] =       gp_weight[p]*J_p*eta[p];

169:     /* form Bt tildeD B */
170:     /*
171:     Ke_ij = Bt_ik . D_kl . B_lj
172:           = B_ki . D_kl . B_lj
173:           = B_ki . D_kk . B_kj
174:     */
175:     for (i = 0; i < 8; i++) {
176:       for (j = 0; j < 8; j++) {
177:         for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
178:           Ke[i+8*j] += B[k][i]*tildeD[k]*B[k][j];
179:         }
180:       }
181:     }
182:   }
183: }

185: static void BForm_Grad(PetscScalar Ke[],PetscScalar coords[])
186: {
187:   PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
188:   PetscScalar Ni_p[NODES_PER_EL],GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
189:   PetscScalar J_p,fac;
190:   PetscInt    p,i,j,di,ngp;

192:   /* define quadrature rule */
193:   CreateGaussQuadrature(&ngp,gp_xi,gp_weight);

195:   /* evaluate bilinear form */
196:   for (p = 0; p < ngp; p++) {
197:     EvaluateBasis_Q1(gp_xi[p],Ni_p);
198:     EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
199:     EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
200:     fac = gp_weight[p]*J_p;

202:     for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
203:       for (di = 0; di < NSD; di++) { /* u dofs */
204:         for (j = 0; j < 4; j++) {  /* p nodes, p dofs = 1 (ie no loop) */
205:           PetscInt IJ;
206:           IJ = map_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);

208:           Ke[IJ] -= GNx_p[di][i]*Ni_p[j]*fac;
209:         }
210:       }
211:     }
212:   }
213: }

215: static void BForm_Div(PetscScalar De[],PetscScalar coords[])
216: {
217:   PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
218:   PetscInt    i,j,nr_g,nc_g;

220:   PetscMemzero(Ge,sizeof(Ge));
221:   BForm_Grad(Ge,coords);

223:   nr_g = U_DOFS*NODES_PER_EL;
224:   nc_g = P_DOFS*NODES_PER_EL;

226:   for (i = 0; i < nr_g; i++) {
227:     for (j = 0; j < nc_g; j++) {
228:       De[nr_g*j+i] = Ge[nc_g*i+j];
229:     }
230:   }
231: }

233: static void BForm_Stabilisation(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
234: {
235:   PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
236:   PetscScalar Ni_p[NODES_PER_EL],GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
237:   PetscScalar J_p,fac,eta_avg;
238:   PetscInt    p,i,j,ngp;

240:   /* define quadrature rule */
241:   CreateGaussQuadrature(&ngp,gp_xi,gp_weight);

243:   /* evaluate bilinear form */
244:   for (p = 0; p < ngp; p++) {
245:     EvaluateBasis_Q1(gp_xi[p],Ni_p);
246:     EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
247:     EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
248:     fac = gp_weight[p]*J_p;

250:     for (i = 0; i < NODES_PER_EL; i++) {
251:       for (j = 0; j < NODES_PER_EL; j++) {
252:         Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.0625);
253:       }
254:     }
255:   }

257:   /* scale */
258:   eta_avg = 0.0;
259:   for (p = 0; p < ngp; p++) eta_avg += eta[p];
260:   eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
261:   fac     = 1.0/eta_avg;
262:   for (i = 0; i < NODES_PER_EL; i++) {
263:     for (j = 0; j < NODES_PER_EL; j++) {
264:       Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
265:     }
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++) {
288:         Ke[NODES_PER_EL*i+j] -= fac*Ni_p[i]*Ni_p[j];
289:       }
290:     }
291:   }

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

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

312:   /* define quadrature rule */
313:   CreateGaussQuadrature(&ngp,gp_xi,gp_weight);

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

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

329: static PetscErrorCode GetElementCoords(const PetscScalar _coords[],const PetscInt e2n[],PetscScalar el_coords[])
330: {
331:   PetscInt i,d;
333:   /* get coords for the element */
334:   for (i=0; i<4; i++) {
335:     for (d=0; d<NSD; d++) {
336:       el_coords[NSD*i+d] = _coords[NSD * e2n[i] + d];
337:     }
338:   }
339:   return(0);
340: }

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

360:   MatZeroEntries(A);
361:   /* setup for coords */
362:   DMGetCoordinateDM(stokes_da,&cda);
363:   DMGetCoordinatesLocal(stokes_da,&coords);
364:   VecGetArrayRead(coords,&_coords);

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

369:   DMDAGetElements(stokes_da,&nel,&npe,&element_list);
370:   for (eidx = 0; eidx < nel; eidx++) {
371:     const PetscInt *element = &element_list[npe*eidx];

373:     /* get coords for the element */
374:     GetElementCoords(_coords,element,el_coords);

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

379:     /* initialise element stiffness matrix */
380:     PetscMemzero(Ae,sizeof(Ae));
381:     PetscMemzero(Ge,sizeof(Ge));
382:     PetscMemzero(De,sizeof(De));
383:     PetscMemzero(Ce,sizeof(Ce));

385:     /* form element stiffness matrix */
386:     BForm_DivT(Ae,el_coords,prop_eta);
387:     BForm_Grad(Ge,el_coords);
388:     BForm_Div(De,el_coords);
389:     BForm_Stabilisation(Ce,el_coords,prop_eta);

391:     /* insert element matrix into global matrix */
392:     DMDAGetElementEqnums_up(element,u_eqn,p_eqn);
393:     MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
394:     MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
395:     MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
396:     MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
397:   }
398:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
399:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

401:   DMSwarmRestoreField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
402:   VecRestoreArrayRead(coords,&_coords);
403:   return(0);
404: }

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

424:   MatZeroEntries(A);
425:   /* setup for coords */
426:   DMGetCoordinateDM(stokes_da,&cda);
427:   DMGetCoordinatesLocal(stokes_da,&coords);
428:   VecGetArrayRead(coords,&_coords);

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

433:   DMDAGetElements(stokes_da,&nel,&npe,&element_list);
434:   for (eidx = 0; eidx < nel; eidx++) {
435:     const PetscInt *element = &element_list[npe*eidx];

437:     /* get coords for the element */
438:     GetElementCoords(_coords,element,el_coords);

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

443:     /* initialise element stiffness matrix */
444:     PetscMemzero(Ae,sizeof(Ae));
445:     PetscMemzero(Ge,sizeof(Ge));
446:     PetscMemzero(De,sizeof(De));
447:     PetscMemzero(Ce,sizeof(Ce));

449:     /* form element stiffness matrix */
450:     BForm_DivT(Ae,el_coords,prop_eta);
451:     BForm_Grad(Ge,el_coords);
452:     BForm_ScaledMassMatrix(Ce,el_coords,prop_eta);

454:     /* insert element matrix into global matrix */
455:     DMDAGetElementEqnums_up(element,u_eqn,p_eqn);
456:     MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
457:     MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
458:     MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
459:     MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
460:   }
461:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
462:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

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

467:   return(0);
468: }

470: static PetscErrorCode AssembleStokes_RHS(Vec F,DM stokes_da,DM quadrature)
471: {
472:   DM                     cda;
473:   Vec                    coords;
474:   const PetscScalar      *_coords;
475:   PetscInt               u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
476:   PetscInt               p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
477:   PetscInt               nel,npe,eidx,i;
478:   const PetscInt         *element_list;
479:   PetscScalar            Fe[NODES_PER_EL*U_DOFS];
480:   PetscScalar            He[NODES_PER_EL*P_DOFS];
481:   PetscScalar            el_coords[NODES_PER_EL*NSD];
482:   PetscScalar            *q_rhs,*prop_fy;
483:   Vec                    local_F;
484:   PetscScalar            *LA_F;
485:   PetscErrorCode         ierr;

488:   VecZeroEntries(F);
489:   /* setup for coords */
490:   DMGetCoordinateDM(stokes_da,&cda);
491:   DMGetCoordinatesLocal(stokes_da,&coords);
492:   VecGetArrayRead(coords,&_coords);

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

497:   /* get acces to the vector */
498:   DMGetLocalVector(stokes_da,&local_F);
499:   VecZeroEntries(local_F);
500:   VecGetArray(local_F,&LA_F);

502:   DMDAGetElements(stokes_da,&nel,&npe,&element_list);
503:   for (eidx = 0; eidx < nel; eidx++) {
504:     const PetscInt *element = &element_list[npe*eidx];

506:     /* get coords for the element */
507:     GetElementCoords(_coords,element,el_coords);

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

512:     /* initialise element stiffness matrix */
513:     PetscMemzero(Fe,sizeof(Fe));
514:     PetscMemzero(He,sizeof(He));

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

519:     /* insert element matrix into global matrix */
520:     DMDAGetElementEqnums_up(element,u_eqn,p_eqn);

522:     for (i=0; i<NODES_PER_EL*U_DOFS; i++) {
523:       LA_F[ u_eqn[i] ] += Fe[i];
524:     }
525:   }
526:   DMSwarmRestoreField(quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
527:   VecRestoreArrayRead(coords,&_coords);

529:   VecRestoreArray(local_F,&LA_F);
530:   DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
531:   DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
532:   DMRestoreLocalVector(stokes_da,&local_F);

534:   return(0);
535: }

537: PetscErrorCode DMSwarmPICInsertPointsCellwise(DM dm,DM dmc,PetscInt e,PetscInt npoints,PetscReal xi[],PetscBool proximity_initialization)
538: {
539:   PetscErrorCode    ierr;
540:   PetscInt          dim,nel,npe,q,k,d,ncurr;
541:   const PetscInt    *element_list;
542:   Vec               coor;
543:   const PetscScalar *_coor;
544:   PetscReal         **basis,*elcoor,*xp;
545:   PetscReal         *swarm_coor;
546:   PetscInt          *swarm_cellid;

549:   DMGetDimension(dm,&dim);
550:   DMDAGetElements(dmc,&nel,&npe,&element_list);

552:   PetscMalloc1(dim*npoints,&xp);
553:   PetscMalloc1(dim*npe,&elcoor);
554:   PetscMalloc1(npoints,&basis);
555:   for (q=0; q<npoints; q++) {
556:     PetscMalloc1(npe,&basis[q]);

558:     switch (dim) {
559:       case 1:
560:         basis[q][0] = 0.5*(1.0 - xi[dim*q+0]);
561:         basis[q][1] = 0.5*(1.0 + xi[dim*q+0]);
562:         break;
563:       case 2:
564:         basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]);
565:         basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]);
566:         basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]);
567:         basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]);
568:         break;

570:       case 3:
571:         basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
572:         basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
573:         basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
574:         basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
575:         basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
576:         basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
577:         basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
578:         basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
579:         break;
580:     }
581:   }

583:   DMGetCoordinatesLocal(dmc,&coor);
584:   VecGetArrayRead(coor,&_coor);
585:   /* compute and store the coordinates for the new points */
586:   {
587:     const PetscInt *element = &element_list[npe*e];

589:     for (k=0; k<npe; k++) {
590:       for (d=0; d<dim; d++) {
591:         elcoor[dim*k+d] = PetscRealPart(_coor[ dim*element[k] + d ]);
592:       }
593:     }
594:     for (q=0; q<npoints; q++) {
595:       for (d=0; d<dim; d++) {
596:         xp[dim*q+d] = 0.0;
597:       }
598:       for (k=0; k<npe; k++) {
599:         for (d=0; d<dim; d++) {
600:           xp[dim*q+d] += basis[q][k] * elcoor[dim*k+d];
601:         }
602:       }
603:     }
604:   }
605:   VecRestoreArrayRead(coor,&_coor);
606:   DMDARestoreElements(dmc,&nel,&npe,&element_list);

608:   DMSwarmGetLocalSize(dm,&ncurr);
609:   DMSwarmAddNPoints(dm,npoints);

611:   if (proximity_initialization) {
612:     PetscInt  *nnlist;
613:     PetscReal *coor_q,*coor_qn;
614:     PetscInt  npoints_e,*plist_e;

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

618:     PetscMalloc1(npoints,&nnlist);
619:     /* find nearest neighour points in this cell */
620:     DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
621:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
622:     for (q=0; q<npoints; q++) {
623:       PetscInt  qn,nearest_neighour = -1;
624:       PetscReal sep,min_sep = PETSC_MAX_REAL;

626:       coor_q = &xp[dim*q];
627:       for (qn=0; qn<npoints_e; qn++) {
628:         coor_qn = &swarm_coor[dim*plist_e[qn]];
629:         sep = 0.0;
630:         for (d=0; d<dim; d++) {
631:           sep += (coor_q[d]-coor_qn[d])*(coor_q[d]-coor_qn[d]);
632:         }
633:         if (sep < min_sep) {
634:           nearest_neighour = plist_e[qn];
635:           min_sep = sep;
636:         }
637:       }
638:       if (nearest_neighour == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cell %D is empty - cannot initalize using nearest neighbours",e);
639:       nnlist[q] = nearest_neighour;
640:     }
641:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
642:     DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

644:     /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
645:     for (q=0; q<npoints; q++) {
646:       DMSwarmCopyPoint(dm,nnlist[q],ncurr+q);
647:     }
648:     DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
649:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
650:     for (q=0; q<npoints; q++) {
651:       /* set the coordinates */
652:       for (d=0; d<dim; d++) {
653:         swarm_coor[dim*(ncurr+q)+d] = xp[dim*q+d];
654:       }
655:       /* set the cell index */
656:       swarm_cellid[ncurr+q] = e;
657:     }
658:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
659:     DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);

661:     PetscFree(plist_e);
662:     PetscFree(nnlist);
663:   } else {
664:     DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
665:     DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
666:     for (q=0; q<npoints; q++) {
667:       /* set the coordinates */
668:       for (d=0; d<dim; d++) {
669:         swarm_coor[dim*(ncurr+q)+d] = xp[dim*q+d];
670:       }
671:       /* set the cell index */
672:       swarm_cellid[ncurr+q] = e;
673:     }
674:     DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
675:     DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
676:   }

678:   PetscFree(xp);
679:   PetscFree(elcoor);
680:   for (q=0; q<npoints; q++) {
681:     PetscFree(basis[q]);
682:   }
683:   PetscFree(basis);
684:   return(0);
685: }

687: PetscErrorCode MaterialPoint_PopulateCell(DM dm_vp,DM dm_mpoint)
688: {
689:   PetscInt        _npe,_nel,e,nel;
690:   const PetscInt  *element;
691:   DM              dmc;
692:   PetscQuadrature quadrature;
693:   const PetscReal *xi;
694:   PetscInt        npoints_q,cnt,cnt_g;
695:   PetscErrorCode  ierr;

698:   DMDAGetElements(dm_vp,&_nel,&_npe,&element);
699:   nel = _nel;
700:   DMDARestoreElements(dm_vp,&_nel,&_npe,&element);

702:   PetscDTGaussTensorQuadrature(2,1,4,-1.0,1.0,&quadrature);
703:   PetscQuadratureGetData(quadrature,NULL,NULL,&npoints_q,&xi,NULL);
704:   DMSwarmGetCellDM(dm_mpoint,&dmc);

706:   DMSwarmSortGetAccess(dm_mpoint);

708:   cnt = 0;
709:   for (e=0; e<nel; e++) {
710:     PetscInt npoints_per_cell;

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

714:     if (npoints_per_cell < 12) {
715:       DMSwarmPICInsertPointsCellwise(dm_mpoint,dm_vp,e,npoints_q,(PetscReal*)xi,PETSC_TRUE);
716:       cnt++;
717:     }
718:   }
719:   MPI_Allreduce(&cnt,&cnt_g,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
720:   if (cnt_g > 0) {
721:     PetscPrintf(PETSC_COMM_WORLD,".... ....pop cont: adjusted %D cells\n",cnt_g);
722:   }

724:   DMSwarmSortRestoreAccess(dm_mpoint);
725:   PetscQuadratureDestroy(&quadrature);
726:   return(0);
727: }

729: PetscErrorCode MaterialPoint_AdvectRK1(DM dm_vp,Vec vp,PetscReal dt,DM dm_mpoint)
730: {
731:   PetscErrorCode    ierr;
732:   Vec               vp_l,coor_l;
733:   const PetscScalar *LA_vp;
734:   PetscInt          i,p,e,npoints,nel,npe;
735:   PetscInt          *mpfield_cell;
736:   PetscReal         *mpfield_coor;
737:   const PetscInt    *element_list;
738:   const PetscInt    *element;
739:   PetscScalar       xi_p[NSD],Ni[NODES_PER_EL];
740:   const PetscScalar *LA_coor;
741:   PetscScalar       dx[NSD];

744:   DMGetCoordinatesLocal(dm_vp,&coor_l);
745:   VecGetArrayRead(coor_l,&LA_coor);

747:   DMGetLocalVector(dm_vp,&vp_l);
748:   DMGlobalToLocalBegin(dm_vp,vp,INSERT_VALUES,vp_l);
749:   DMGlobalToLocalEnd(dm_vp,vp,INSERT_VALUES,vp_l);
750:   VecGetArrayRead(vp_l,&LA_vp);

752:   DMDAGetElements(dm_vp,&nel,&npe,&element_list);
753:   DMSwarmGetLocalSize(dm_mpoint,&npoints);
754:   DMSwarmGetField(dm_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
755:   DMSwarmGetField(dm_mpoint,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
756:   for (p=0; p<npoints; p++) {
757:     PetscReal         *coor_p;
758:     PetscScalar       vel_n[NSD*NODES_PER_EL],vel_p[NSD];
759:     const PetscScalar *x0;
760:     const PetscScalar *x2;

762:     e       = mpfield_cell[p];
763:     coor_p  = &mpfield_coor[NSD*p];
764:     element = &element_list[NODES_PER_EL*e];

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

770:     dx[0] = x2[0] - x0[0];
771:     dx[1] = x2[1] - x0[1];

773:     xi_p[0] = 2.0 * (coor_p[0] - x0[0])/dx[0] - 1.0;
774:     xi_p[1] = 2.0 * (coor_p[1] - x0[1])/dx[1] - 1.0;
775:     if (PetscRealPart(xi_p[0]) < -1.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"value (xi) too small %1.4e [e=%D]\n",(double)PetscRealPart(xi_p[0]),e);
776:     if (PetscRealPart(xi_p[0]) >  1.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"value (xi) too large %1.4e [e=%D]\n",(double)PetscRealPart(xi_p[0]),e);
777:     if (PetscRealPart(xi_p[1]) < -1.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"value (eta) too small %1.4e [e=%D]\n",(double)PetscRealPart(xi_p[1]),e);
778:     if (PetscRealPart(xi_p[1]) >  1.0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"value (eta) too large %1.4e [e=%D]\n",(double)PetscRealPart(xi_p[1]),e);

780:     /* evaluate basis functions */
781:     EvaluateBasis_Q1(xi_p,Ni);

783:     /* get cell nodal velocities */
784:     for (i=0; i<NODES_PER_EL; i++) {
785:       PetscInt nid;

787:       nid = element[i];
788:       vel_n[NSD*i+0] = LA_vp[(NSD+1)*nid+0];
789:       vel_n[NSD*i+1] = LA_vp[(NSD+1)*nid+1];
790:     }

792:     /* interpolate velocity */
793:     vel_p[0] = vel_p[1] = 0.0;
794:     for (i=0; i<NODES_PER_EL; i++) {
795:       vel_p[0] += Ni[i] * vel_n[NSD*i+0];
796:       vel_p[1] += Ni[i] * vel_n[NSD*i+1];
797:     }

799:     coor_p[0] += dt * PetscRealPart(vel_p[0]);
800:     coor_p[1] += dt * PetscRealPart(vel_p[1]);
801:   }

803:   DMSwarmRestoreField(dm_mpoint,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
804:   DMSwarmRestoreField(dm_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
805:   DMDARestoreElements(dm_vp,&nel,&npe,&element_list);
806:   VecRestoreArrayRead(vp_l,&LA_vp);
807:   DMRestoreLocalVector(dm_vp,&vp_l);
808:   VecRestoreArrayRead(coor_l,&LA_coor);
809:   return(0);
810: }

812: PetscErrorCode MaterialPoint_Interpolate(DM dm,Vec eta_v,Vec rho_v,DM dm_quadrature)
813: {
814:   Vec            eta_l,rho_l;
815:   PetscScalar    *_eta_l,*_rho_l;
816:   PetscInt       nqp,npe,nel;
817:   PetscScalar    qp_xi[GAUSS_POINTS][NSD];
818:   PetscScalar    qp_weight[GAUSS_POINTS];
819:   PetscInt       q,k,e;
820:   PetscScalar    Ni[GAUSS_POINTS][NODES_PER_EL];
821:   const PetscInt *element_list;
822:   PetscReal      *q_eta,*q_rhs;

826:   /* define quadrature rule */
827:   CreateGaussQuadrature(&nqp,qp_xi,qp_weight);
828:   for (q=0; q<nqp; q++) {
829:     EvaluateBasis_Q1(qp_xi[q],Ni[q]);
830:   }

832:   DMGetLocalVector(dm,&eta_l);
833:   DMGetLocalVector(dm,&rho_l);

835:   DMGlobalToLocalBegin(dm,eta_v,INSERT_VALUES,eta_l);
836:   DMGlobalToLocalEnd(dm,eta_v,INSERT_VALUES,eta_l);
837:   DMGlobalToLocalBegin(dm,rho_v,INSERT_VALUES,rho_l);
838:   DMGlobalToLocalEnd(dm,rho_v,INSERT_VALUES,rho_l);

840:   VecGetArray(eta_l,&_eta_l);
841:   VecGetArray(rho_l,&_rho_l);

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

846:   DMDAGetElements(dm,&nel,&npe,&element_list);
847:   for (e=0; e<nel; e++) {
848:     PetscScalar    eta_field_e[NODES_PER_EL];
849:     PetscScalar    rho_field_e[NODES_PER_EL];
850:     const PetscInt *element = &element_list[4*e];

852:     for (k=0; k<NODES_PER_EL; k++) {
853:       eta_field_e[k] = _eta_l[ element[k] ];
854:       rho_field_e[k] = _rho_l[ element[k] ];
855:     }

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

860:       eta_q = rho_q = 0.0;
861:       for (k=0; k<NODES_PER_EL; k++) {
862:         eta_q += Ni[q][k] * eta_field_e[k];
863:         rho_q += Ni[q][k] * rho_field_e[k];
864:       }

866:       q_eta[nqp*e+q] = PetscRealPart(eta_q);
867:       q_rhs[nqp*e+q] = PetscRealPart(rho_q);
868:     }
869:   }
870:   DMDARestoreElements(dm,&nel,&npe,&element_list);

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

875:   VecRestoreArray(rho_l,&_rho_l);
876:   VecRestoreArray(eta_l,&_eta_l);
877:   DMRestoreLocalVector(dm,&rho_l);
878:   DMRestoreLocalVector(dm,&eta_l);
879:   return(0);
880: }

882: static PetscErrorCode SolveTimeDepStokes(PetscInt mx,PetscInt my)
883: {
884:   DM                     dm_stokes,dm_coeff;
885:   PetscInt               u_dof,p_dof,dof,stencil_width;
886:   Mat                    A,B;
887:   PetscInt               nel_local;
888:   Vec                    eta_v,rho_v;
889:   Vec                    f,X;
890:   KSP                    ksp;
891:   PC                     pc;
892:   char                   filename[PETSC_MAX_PATH_LEN];
893:   DM                     dms_quadrature,dms_mpoint;
894:   PetscInt               nel,npe,npoints;
895:   const PetscInt         *element_list;
896:   PetscInt               tk,nt,dump_freq;
897:   PetscReal              dt,dt_max = 0.0;
898:   PetscReal              vx[2],vy[2],max_v = 0.0,max_v_step,dh;
899:   PetscErrorCode         ierr;
900:   const char             *fieldnames[] = { "eta" , "rho" };
901:   Vec                    *pfields;
902:   PetscInt               ppcell = 1;
903:   PetscReal              time,delta_eta = 1.0;
904:   PetscBool              randomize_coords = PETSC_FALSE;
905:   PetscReal              randomize_fac = 0.25;
906:   PetscBool              no_view = PETSC_FALSE;
907:   PetscBool              isbddc;

910:   /*
911:     Generate the DMDA for the velocity and pressure spaces.
912:     We use Q1 elements for both fields.
913:     The Q1 FE basis on a regular mesh has a 9-point stencil (DMDA_STENCIL_BOX)
914:     The number of nodes in each direction is mx+1, my+1
915:   */
916:   u_dof         = U_DOFS; /* Vx, Vy - velocities */
917:   p_dof         = P_DOFS; /* p - pressure */
918:   dof           = u_dof + p_dof;
919:   stencil_width = 1;
920:   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);
921:   DMDASetElementType(dm_stokes,DMDA_ELEMENT_Q1);
922:   DMSetMatType(dm_stokes,MATAIJ);
923:   DMSetFromOptions(dm_stokes);
924:   DMSetUp(dm_stokes);
925:   DMDASetFieldName(dm_stokes,0,"ux");
926:   DMDASetFieldName(dm_stokes,1,"uy");
927:   DMDASetFieldName(dm_stokes,2,"p");

929:   /* unit box [0,0.9142] x [0,1] */
930:   DMDASetUniformCoordinates(dm_stokes,0.0,0.9142,0.0,1.0,0.,0.);
931:   dh = 1.0/((PetscReal)(mx));

933:   /* Get local number of elements */
934:   {
935:     DMDAGetElements(dm_stokes,&nel,&npe,&element_list);

937:     nel_local = nel;

939:     DMDARestoreElements(dm_stokes,&nel,&npe,&element_list);
940:   }

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

945:   /* Create the swarm for storing quadrature point values */
946:   DMCreate(PETSC_COMM_WORLD,&dms_quadrature);
947:   DMSetType(dms_quadrature,DMSWARM);
948:   DMSetDimension(dms_quadrature,2);

950:   /* Register fields for viscosity and density on the quadrature points */
951:   DMSwarmRegisterPetscDatatypeField(dms_quadrature,"eta_q",1,PETSC_REAL);
952:   DMSwarmRegisterPetscDatatypeField(dms_quadrature,"rho_q",1,PETSC_REAL);
953:   DMSwarmFinalizeFieldRegister(dms_quadrature);
954:   DMSwarmSetLocalSizes(dms_quadrature,nel_local * GAUSS_POINTS,0);

956:   /* Create the material point swarm */
957:   DMCreate(PETSC_COMM_WORLD,&dms_mpoint);
958:   DMSetType(dms_mpoint,DMSWARM);
959:   DMSetDimension(dms_mpoint,2);

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

964:   /*
965:      Specify the DM to use for point location and projections
966:      within the context of a PIC scheme
967:   */
968:   DMSwarmSetCellDM(dms_mpoint,dm_coeff);

970:   /* Register fields for viscosity and density */
971:   DMSwarmRegisterPetscDatatypeField(dms_mpoint,"eta",1,PETSC_REAL);
972:   DMSwarmRegisterPetscDatatypeField(dms_mpoint,"rho",1,PETSC_REAL);
973:   DMSwarmFinalizeFieldRegister(dms_mpoint);

975:   PetscOptionsGetInt(NULL,NULL,"-ppcell",&ppcell,NULL);
976:   DMSwarmSetLocalSizes(dms_mpoint,nel_local * ppcell,100);

978:   /*
979:     Layout the material points in space using the cell DM.
980:     Particle coordinates are defined by cell wise using different methods.
981:     - DMSWARMPIC_LAYOUT_GAUSS defines particles coordinates at the positions
982:                               corresponding to a Gauss quadrature rule with
983:                               ppcell points in each direction.
984:     - DMSWARMPIC_LAYOUT_REGULAR defines particle coordinates at the centoid of
985:                                 ppcell x ppcell quadralaterals defined within the
986:                                 reference element.
987:     - DMSWARMPIC_LAYOUT_SUBDIVISION defines particles coordinates at the centroid
988:                                     of each quadralateral obtained by sub-dividing
989:                                     the reference element cell ppcell times.
990:   */
991:   DMSwarmInsertPointsUsingCellDM(dms_mpoint,DMSWARMPIC_LAYOUT_SUBDIVISION,ppcell);

993:   /*
994:     Defne a high resolution layer of material points across the material interface
995:   */
996:   {
997:     PetscInt  npoints_dir_x[2];
998:     PetscReal min[2],max[2];

1000:     npoints_dir_x[0] = (PetscInt)(0.9142/(0.05*dh));
1001:     npoints_dir_x[1] = (PetscInt)((0.25-0.15)/(0.05*dh));
1002:     min[0] = 0.0;  max[0] = 0.9142;
1003:     min[1] = 0.05; max[1] = 0.35;
1004:     DMSwarmSetPointsUniformCoordinates(dms_mpoint,min,max,npoints_dir_x,ADD_VALUES);
1005:   }

1007:   /*
1008:     Define a high resolution layer of material points near the surface of the domain
1009:     to deal with weakly compressible Q1-Q1 elements. These elements "self compact"
1010:     when applied to buouyancy driven flow. The error in div(u) is O(h).
1011:   */
1012:   {
1013:     PetscInt  npoints_dir_x[2];
1014:     PetscReal min[2],max[2];

1016:     npoints_dir_x[0] = (PetscInt)(0.9142/(0.25*dh));
1017:     npoints_dir_x[1] = (PetscInt)(3.0*dh/(0.25*dh));
1018:     min[0] = 0.0;          max[0] = 0.9142;
1019:     min[1] = 1.0 - 3.0*dh; max[1] = 1.0-0.0001;
1020:     DMSwarmSetPointsUniformCoordinates(dms_mpoint,min,max,npoints_dir_x,ADD_VALUES);
1021:   }

1023:   DMView(dms_mpoint,PETSC_VIEWER_STDOUT_WORLD);

1025:   /* Define initial material properties on each particle in the material point swarm */
1026:   PetscOptionsGetReal(NULL,NULL,"-delta_eta",&delta_eta,NULL);
1027:   PetscOptionsGetBool(NULL,NULL,"-randomize_coords",&randomize_coords,NULL);
1028:   PetscOptionsGetReal(NULL,NULL,"-randomize_fac",&randomize_fac,NULL);
1029:   if (randomize_fac > 1.0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"The value of -randomize_fac should be <= 1.0");
1030:   {
1031:     PetscReal   *array_x,*array_e,*array_r;
1032:     PetscInt    p;
1033:     PetscRandom r;
1034:     PetscMPIInt rank;

1036:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

1038:     PetscRandomCreate(PETSC_COMM_SELF,&r);
1039:     PetscRandomSetInterval(r,-randomize_fac*dh,randomize_fac*dh);
1040:     PetscRandomSetSeed(r,(unsigned long)rank);
1041:     PetscRandomSeed(r);

1043:     DMDAGetElements(dm_stokes,&nel,&npe,&element_list);

1045:     /*
1046:        Fetch the registered data from the material point DMSwarm.
1047:        The fields "eta" and "rho" were registered by this example.
1048:        The field identified by the the variable DMSwarmPICField_coor
1049:        was registered by the DMSwarm implementation when the function
1050:          DMSwarmSetType(dms_mpoint,DMSWARM_PIC)
1051:        was called. The returned array defines the coordinates of each
1052:        material point in the point swarm.
1053:     */
1054:     DMSwarmGetField(dms_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&array_x);
1055:     DMSwarmGetField(dms_mpoint,"eta",               NULL,NULL,(void**)&array_e);
1056:     DMSwarmGetField(dms_mpoint,"rho",               NULL,NULL,(void**)&array_r);

1058:     DMSwarmGetLocalSize(dms_mpoint,&npoints);
1059:     for (p = 0; p < npoints; p++) {
1060:       PetscReal x_p[2],rr[2];

1062:       if (randomize_coords) {
1063:         PetscRandomGetValueReal(r,&rr[0]);
1064:         PetscRandomGetValueReal(r,&rr[1]);
1065:         array_x[2*p + 0] += rr[0];
1066:         array_x[2*p + 1] += rr[1];
1067:       }

1069:       /* Get the coordinates of point, p */
1070:       x_p[0] = array_x[2*p + 0];
1071:       x_p[1] = array_x[2*p + 1];

1073:        if (x_p[1] < (0.2 + 0.02*PetscCosReal(PETSC_PI*x_p[0]/0.9142))) {
1074:          /* Material properties below the interface */
1075:          array_e[p] = 1.0 * (1.0/delta_eta);
1076:          array_r[p] = 0.0;
1077:        } else {
1078:          /* Material properties above the interface */
1079:          array_e[p] = 1.0;
1080:          array_r[p] = 1.0;
1081:        }
1082:     }

1084:     /*
1085:        Restore the fetched data fields from the material point DMSwarm.
1086:        Calling the Restore function invalidates the points array_r, array_e, array_x
1087:        by setting them to NULL.
1088:     */
1089:     DMSwarmRestoreField(dms_mpoint,"rho",NULL,NULL,(void**)&array_r);
1090:     DMSwarmRestoreField(dms_mpoint,"eta",NULL,NULL,(void**)&array_e);
1091:     DMSwarmRestoreField(dms_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&array_x);

1093:     DMDARestoreElements(dm_stokes,&nel,&npe,&element_list);
1094:     PetscRandomDestroy(&r);
1095:   }

1097:   /*
1098:      If the particle coordinates where randomly shifted, they may have crossed into another
1099:      element, or into another sub-domain. To account for this we call the Migrate function.
1100:   */
1101:   if (randomize_coords) {
1102:     DMSwarmMigrate(dms_mpoint,PETSC_TRUE);
1103:   }

1105:   PetscOptionsGetBool(NULL,NULL,"-no_view",&no_view,NULL);
1106:   if (!no_view) {
1107:     DMSwarmViewXDMF(dms_mpoint,"ic_coeff_dms.xmf");
1108:   }

1110:   /* project the swarm properties */
1111:   DMSwarmProjectFields(dms_mpoint,2,fieldnames,&pfields,PETSC_FALSE);
1112:   eta_v = pfields[0];
1113:   rho_v = pfields[1];
1114:   PetscObjectSetName((PetscObject)eta_v,"eta");
1115:   PetscObjectSetName((PetscObject)rho_v,"rho");
1116:   MaterialPoint_Interpolate(dm_coeff,eta_v,rho_v,dms_quadrature);

1118:   /* view projected coefficients eta and rho */
1119:   if (!no_view) {
1120:     PetscViewer viewer;

1122:     PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1123:     PetscViewerSetType(viewer,PETSCVIEWERVTK);
1124:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1125:     PetscViewerFileSetName(viewer,"ic_coeff_dmda.vts");
1126:     VecView(eta_v,viewer);
1127:     VecView(rho_v,viewer);
1128:     PetscViewerDestroy(&viewer);
1129:   }

1131:   DMCreateMatrix(dm_stokes,&A);
1132:   DMCreateMatrix(dm_stokes,&B);
1133:   DMCreateGlobalVector(dm_stokes,&f);
1134:   DMCreateGlobalVector(dm_stokes,&X);

1136:   AssembleStokes_A(A,dm_stokes,dms_quadrature);
1137:   AssembleStokes_PC(B,dm_stokes,dms_quadrature);
1138:   AssembleStokes_RHS(f,dm_stokes,dms_quadrature);

1140:   DMDAApplyBoundaryConditions(dm_stokes,A,f);
1141:   DMDAApplyBoundaryConditions(dm_stokes,B,NULL);

1143:   KSPCreate(PETSC_COMM_WORLD,&ksp);
1144:   KSPSetOptionsPrefix(ksp,"stokes_");
1145:   KSPSetDM(ksp,dm_stokes);
1146:   KSPSetDMActive(ksp,PETSC_FALSE);
1147:   KSPSetOperators(ksp,A,B);
1148:   KSPSetFromOptions(ksp);
1149:   KSPGetPC(ksp,&pc);
1150:   PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);
1151:   if (isbddc) {
1152:     KSPSetOperators(ksp,A,A);
1153:   }

1155:   /* Define u-v-p indices for fieldsplit */
1156:   {
1157:     PC             pc;
1158:     const PetscInt ufields[] = {0,1},pfields[1] = {2};

1160:     KSPGetPC(ksp,&pc);
1161:     PCFieldSplitSetBlockSize(pc,3);
1162:     PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1163:     PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1164:   }

1166:   /* If using a fieldsplit preconditioner, attach a DMDA to the velocity split so that geometric multigrid can be used */
1167:   {
1168:     PC        pc,pc_u;
1169:     KSP       *sub_ksp,ksp_u;
1170:     PetscInt  nsplits;
1171:     DM        dm_u;
1172:     PetscBool is_pcfs;

1174:     KSPGetPC(ksp,&pc);

1176:     is_pcfs = PETSC_FALSE;
1177:     PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&is_pcfs);

1179:     if (is_pcfs) {
1180:       KSPSetUp(ksp);
1181:       KSPGetPC(ksp,&pc);
1182:       PCFieldSplitGetSubKSP(pc,&nsplits,&sub_ksp);
1183:       ksp_u = sub_ksp[0];
1184:       PetscFree(sub_ksp);

1186:       if (nsplits == 2) {
1187:         DMDACreateCompatibleDMDA(dm_stokes,2,&dm_u);

1189:         KSPSetDM(ksp_u,dm_u);
1190:         KSPSetDMActive(ksp_u,PETSC_FALSE);
1191:         DMDestroy(&dm_u);

1193:         /* enforce galerkin coarse grids be used */
1194:         KSPGetPC(ksp_u,&pc_u);
1195:         PCMGSetGalerkin(pc_u,PC_MG_GALERKIN_PMAT);
1196:       }
1197:     }
1198:   }

1200:   dump_freq = 10;
1201:   PetscOptionsGetInt(NULL,NULL,"-dump_freq",&dump_freq,NULL);
1202:   nt = 10;
1203:   PetscOptionsGetInt(NULL,NULL,"-nt",&nt,NULL);
1204:   time = 0.0;
1205:   for (tk=1; tk<=nt; tk++) {

1207:     PetscPrintf(PETSC_COMM_WORLD,".... assemble\n");
1208:     AssembleStokes_A(A,dm_stokes,dms_quadrature);
1209:     AssembleStokes_PC(B,dm_stokes,dms_quadrature);
1210:     AssembleStokes_RHS(f,dm_stokes,dms_quadrature);

1212:     PetscPrintf(PETSC_COMM_WORLD,".... bc imposition\n");
1213:     DMDAApplyBoundaryConditions(dm_stokes,A,f);
1214:     DMDAApplyBoundaryConditions(dm_stokes,B,NULL);

1216:     PetscPrintf(PETSC_COMM_WORLD,".... solve\n");
1217:     KSPSetOperators(ksp,A, isbddc ? A : B);
1218:     KSPSolve(ksp,f,X);

1220:     VecStrideMax(X,0,NULL,&vx[1]);
1221:     VecStrideMax(X,1,NULL,&vy[1]);
1222:     VecStrideMin(X,0,NULL,&vx[0]);
1223:     VecStrideMin(X,1,NULL,&vy[0]);

1225:     max_v_step = PetscMax(vx[0],vx[1]);
1226:     max_v_step = PetscMax(max_v_step,vy[0]);
1227:     max_v_step = PetscMax(max_v_step,vy[1]);
1228:     max_v = PetscMax(max_v,max_v_step);

1230:     dt_max = 2.0;
1231:     dt = 0.5 * (dh / max_v_step);
1232:     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);
1233:     dt = PetscMin(dt_max,dt);

1235:     /* advect */
1236:     PetscPrintf(PETSC_COMM_WORLD,".... advect\n");
1237:     MaterialPoint_AdvectRK1(dm_stokes,X,dt,dms_mpoint);

1239:     /* migrate */
1240:     PetscPrintf(PETSC_COMM_WORLD,".... migrate\n");
1241:     DMSwarmMigrate(dms_mpoint,PETSC_TRUE);

1243:     /* update cell population */
1244:     PetscPrintf(PETSC_COMM_WORLD,".... populate cells\n");
1245:     MaterialPoint_PopulateCell(dm_stokes,dms_mpoint);

1247:     /* update coefficients on quadrature points */
1248:     PetscPrintf(PETSC_COMM_WORLD,".... project\n");
1249:     DMSwarmProjectFields(dms_mpoint,2,fieldnames,&pfields,PETSC_TRUE);
1250:     eta_v = pfields[0];
1251:     rho_v = pfields[1];
1252:     PetscPrintf(PETSC_COMM_WORLD,".... interp\n");
1253:     MaterialPoint_Interpolate(dm_coeff,eta_v,rho_v,dms_quadrature);

1255:     if (tk%dump_freq == 0) {
1256:       PetscViewer viewer;

1258:       PetscPrintf(PETSC_COMM_WORLD,".... write XDMF, VTS\n");
1259:       PetscSNPrintf(filename,PETSC_MAX_PATH_LEN-1,"step%.4D_coeff_dms.xmf",tk);
1260:       DMSwarmViewXDMF(dms_mpoint,filename);

1262:       PetscSNPrintf(filename,PETSC_MAX_PATH_LEN-1,"step%.4D_vp_dm.vts",tk);
1263:       PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1264:       PetscViewerSetType(viewer,PETSCVIEWERVTK);
1265:       PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1266:       PetscViewerFileSetName(viewer,filename);
1267:       VecView(X,viewer);
1268:       PetscViewerDestroy(&viewer);
1269:     }
1270:     time += dt;
1271:     PetscPrintf(PETSC_COMM_WORLD,"step %D : time %1.2e \n",tk,(double)time);
1272:   }

1274:   KSPDestroy(&ksp);
1275:   VecDestroy(&X);
1276:   VecDestroy(&f);
1277:   MatDestroy(&A);
1278:   MatDestroy(&B);
1279:   VecDestroy(&eta_v);
1280:   VecDestroy(&rho_v);
1281:   PetscFree(pfields);

1283:   DMDestroy(&dms_mpoint);
1284:   DMDestroy(&dms_quadrature);
1285:   DMDestroy(&dm_coeff);
1286:   DMDestroy(&dm_stokes);
1287:   return(0);
1288: }

1290: /*
1291:  <sequential run>
1292:  ./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
1293: */
1294: int main(int argc,char **args)
1295: {
1297:   PetscInt       mx,my;
1298:   PetscBool      set = PETSC_FALSE;

1300:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
1301:   mx = my = 10;
1302:   PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1303:   PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1304:   PetscOptionsGetInt(NULL,NULL,"-mxy",&mx,&set);
1305:   if (set) {
1306:     my = mx;
1307:   }
1308:   SolveTimeDepStokes(mx,my);
1309:   PetscFinalize();
1310:   return ierr;
1311: }

1313: /* -------------------------- helpers for boundary conditions -------------------------------- */
1314: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
1315: {
1316:   DM                     cda;
1317:   Vec                    coords;
1318:   PetscInt               si,sj,nx,ny,i,j;
1319:   PetscInt               M,N;
1320:   DMDACoor2d             **_coords;
1321:   const PetscInt         *g_idx;
1322:   PetscInt               *bc_global_ids;
1323:   PetscScalar            *bc_vals;
1324:   PetscInt               nbcs;
1325:   PetscInt               n_dofs;
1326:   PetscErrorCode         ierr;
1327:   ISLocalToGlobalMapping ltogm;

1330:   DMGetLocalToGlobalMapping(da,&ltogm);
1331:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1333:   DMGetCoordinateDM(da,&cda);
1334:   DMGetCoordinatesLocal(da,&coords);
1335:   DMDAVecGetArray(cda,coords,&_coords);
1336:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1337:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1339:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
1340:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

1345:   i = nx-1;
1346:   for (j = 0; j < ny; j++) {
1347:     PetscInt local_id;

1349:     local_id = i+j*nx;

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

1353:     bc_vals[j] =  0.0;
1354:   }
1355:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1356:   nbcs = 0;
1357:   if ((si+nx) == (M)) nbcs = ny;

1359:   if (b) {
1360:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1361:     VecAssemblyBegin(b);
1362:     VecAssemblyEnd(b);
1363:   }
1364:   if (A) {
1365:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1366:   }

1368:   PetscFree(bc_vals);
1369:   PetscFree(bc_global_ids);

1371:   DMDAVecRestoreArray(cda,coords,&_coords);
1372:   return(0);
1373: }

1375: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
1376: {
1377:   DM                     cda;
1378:   Vec                    coords;
1379:   PetscInt               si,sj,nx,ny,i,j;
1380:   PetscInt               M,N;
1381:   DMDACoor2d             **_coords;
1382:   const PetscInt         *g_idx;
1383:   PetscInt               *bc_global_ids;
1384:   PetscScalar            *bc_vals;
1385:   PetscInt               nbcs;
1386:   PetscInt               n_dofs;
1387:   PetscErrorCode         ierr;
1388:   ISLocalToGlobalMapping ltogm;

1391:   DMGetLocalToGlobalMapping(da,&ltogm);
1392:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1394:   DMGetCoordinateDM(da,&cda);
1395:   DMGetCoordinatesLocal(da,&coords);
1396:   DMDAVecGetArray(cda,coords,&_coords);
1397:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1398:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1400:   PetscMalloc1(ny*n_dofs,&bc_global_ids);
1401:   PetscMalloc1(ny*n_dofs,&bc_vals);

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

1406:   i = 0;
1407:   for (j = 0; j < ny; j++) {
1408:     PetscInt local_id;

1410:     local_id = i+j*nx;

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

1414:     bc_vals[j] =  0.0;
1415:   }
1416:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1417:   nbcs = 0;
1418:   if (si == 0) nbcs = ny;

1420:   if (b) {
1421:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1422:     VecAssemblyBegin(b);
1423:     VecAssemblyEnd(b);
1424:   }

1426:   if (A) {
1427:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1428:   }

1430:   PetscFree(bc_vals);
1431:   PetscFree(bc_global_ids);

1433:   DMDAVecRestoreArray(cda,coords,&_coords);
1434:   return(0);
1435: }

1437: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
1438: {
1439:   DM                     cda;
1440:   Vec                    coords;
1441:   PetscInt               si,sj,nx,ny,i,j;
1442:   PetscInt               M,N;
1443:   DMDACoor2d             **_coords;
1444:   const PetscInt         *g_idx;
1445:   PetscInt               *bc_global_ids;
1446:   PetscScalar            *bc_vals;
1447:   PetscInt               nbcs;
1448:   PetscInt               n_dofs;
1449:   PetscErrorCode         ierr;
1450:   ISLocalToGlobalMapping ltogm;

1453:   DMGetLocalToGlobalMapping(da,&ltogm);
1454:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1456:   DMGetCoordinateDM(da,&cda);
1457:   DMGetCoordinatesLocal(da,&coords);
1458:   DMDAVecGetArray(cda,coords,&_coords);
1459:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1460:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1462:   PetscMalloc1(nx,&bc_global_ids);
1463:   PetscMalloc1(nx,&bc_vals);

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

1468:   j = ny-1;
1469:   for (i = 0; i < nx; i++) {
1470:     PetscInt local_id;

1472:     local_id = i+j*nx;

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

1476:     bc_vals[i] =  0.0;
1477:   }
1478:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1479:   nbcs = 0;
1480:   if ((sj+ny) == (N)) nbcs = nx;

1482:   if (b) {
1483:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1484:     VecAssemblyBegin(b);
1485:     VecAssemblyEnd(b);
1486:   }
1487:   if (A) {
1488:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,NULL,NULL);
1489:   }

1491:   PetscFree(bc_vals);
1492:   PetscFree(bc_global_ids);

1494:   DMDAVecRestoreArray(cda,coords,&_coords);
1495:   return(0);
1496: }

1498: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
1499: {
1500:   DM                     cda;
1501:   Vec                    coords;
1502:   PetscInt               si,sj,nx,ny,i,j;
1503:   PetscInt               M,N;
1504:   DMDACoor2d             **_coords;
1505:   const PetscInt         *g_idx;
1506:   PetscInt               *bc_global_ids;
1507:   PetscScalar            *bc_vals;
1508:   PetscInt               nbcs;
1509:   PetscInt               n_dofs;
1510:   PetscErrorCode         ierr;
1511:   ISLocalToGlobalMapping ltogm;

1514:   DMGetLocalToGlobalMapping(da,&ltogm);
1515:   ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);

1517:   DMGetCoordinateDM(da,&cda);
1518:   DMGetCoordinatesLocal(da,&coords);
1519:   DMDAVecGetArray(cda,coords,&_coords);
1520:   DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1521:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);

1523:   PetscMalloc1(nx,&bc_global_ids);
1524:   PetscMalloc1(nx,&bc_vals);

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

1529:   j = 0;
1530:   for (i = 0; i < nx; i++) {
1531:     PetscInt local_id;

1533:     local_id = i+j*nx;

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

1537:     bc_vals[i] =  0.0;
1538:   }
1539:   ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1540:   nbcs = 0;
1541:   if (sj == 0) nbcs = nx;

1543:   if (b) {
1544:     VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1545:     VecAssemblyBegin(b);
1546:     VecAssemblyEnd(b);
1547:   }
1548:   if (A) {
1549:     MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1550:   }

1552:   PetscFree(bc_vals);
1553:   PetscFree(bc_global_ids);

1555:   DMDAVecRestoreArray(cda,coords,&_coords);
1556:   return(0);
1557: }

1559: /*
1560:  Impose free slip boundary conditions on the left/right faces: u_i n_i = 0, tau_{ij} t_j = 0
1561:  Impose no slip boundray conditions on the top/bottom faces:   u_i n_i = 0, u_i t_i = 0
1562: */
1563: static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes,Mat A,Vec f)
1564: {

1568:   BCApplyZero_NORTH(dm_stokes,0,A,f);
1569:   BCApplyZero_NORTH(dm_stokes,1,A,f);
1570:   BCApplyZero_EAST(dm_stokes,0,A,f);
1571:   BCApplyZero_SOUTH(dm_stokes,0,A,f);
1572:   BCApplyZero_SOUTH(dm_stokes,1,A,f);
1573:   BCApplyZero_WEST(dm_stokes,0,A,f);
1574:   return(0);
1575: }

1577: /*TEST

1579:    test:
1580:      suffix: 1
1581:      args: -no_view
1582:      requires: !complex double
1583:      filter: grep -v atomic
1584:      filter_output: grep -v atomic
1585:    test:
1586:      suffix: 1_matis
1587:      requires: !complex double
1588:      args: -no_view -dm_mat_type is
1589:      filter: grep -v atomic
1590:      filter_output: grep -v atomic
1591:    testset:
1592:      nsize: 4
1593:      requires: !complex double
1594:      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
1595:      filter: grep -v atomic
1596:      filter_output: grep -v atomic
1597:      test:
1598:        suffix: fetidp
1599:        args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
1600:      test:
1601:        suffix: fetidp_lumped
1602:        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
1603:      test:
1604:        suffix: fetidp_saddlepoint
1605:        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
1606:      test:
1607:        suffix: fetidp_saddlepoint_lumped
1608:        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
1609: TEST*/