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;
359: MatZeroEntries(A);
360: /* setup for coords */
361: DMGetCoordinateDM(stokes_da,&cda);
362: DMGetCoordinatesLocal(stokes_da,&coords);
363: VecGetArrayRead(coords,&_coords);
365: /* setup for coefficients */
366: DMSwarmGetField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
368: DMDAGetElements(stokes_da,&nel,&npe,&element_list);
369: for (eidx = 0; eidx < nel; eidx++) {
370: const PetscInt *element = &element_list[npe*eidx];
372: /* get coords for the element */
373: GetElementCoords(_coords,element,el_coords);
375: /* get coefficients for the element */
376: prop_eta = &q_eta[GAUSS_POINTS * eidx];
378: /* initialise element stiffness matrix */
379: PetscMemzero(Ae,sizeof(Ae));
380: PetscMemzero(Ge,sizeof(Ge));
381: PetscMemzero(De,sizeof(De));
382: PetscMemzero(Ce,sizeof(Ce));
384: /* form element stiffness matrix */
385: BForm_DivT(Ae,el_coords,prop_eta);
386: BForm_Grad(Ge,el_coords);
387: BForm_Div(De,el_coords);
388: BForm_Stabilisation(Ce,el_coords,prop_eta);
390: /* insert element matrix into global matrix */
391: DMDAGetElementEqnums_up(element,u_eqn,p_eqn);
392: MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
393: MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
394: MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
395: MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
396: }
397: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
398: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
400: DMSwarmRestoreField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
401: VecRestoreArrayRead(coords,&_coords);
402: return 0;
403: }
405: static PetscErrorCode AssembleStokes_PC(Mat A,DM stokes_da,DM quadrature)
406: {
407: DM cda;
408: Vec coords;
409: const PetscScalar *_coords;
410: PetscInt u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
411: PetscInt p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
412: PetscInt nel,npe,eidx;
413: const PetscInt *element_list;
414: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
415: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
416: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
417: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
418: PetscScalar el_coords[NODES_PER_EL*NSD];
419: PetscScalar *q_eta,*prop_eta;
422: MatZeroEntries(A);
423: /* setup for coords */
424: DMGetCoordinateDM(stokes_da,&cda);
425: DMGetCoordinatesLocal(stokes_da,&coords);
426: VecGetArrayRead(coords,&_coords);
428: /* setup for coefficients */
429: DMSwarmGetField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
431: DMDAGetElements(stokes_da,&nel,&npe,&element_list);
432: for (eidx = 0; eidx < nel; eidx++) {
433: const PetscInt *element = &element_list[npe*eidx];
435: /* get coords for the element */
436: GetElementCoords(_coords,element,el_coords);
438: /* get coefficients for the element */
439: prop_eta = &q_eta[GAUSS_POINTS * eidx];
441: /* initialise element stiffness matrix */
442: PetscMemzero(Ae,sizeof(Ae));
443: PetscMemzero(Ge,sizeof(Ge));
444: PetscMemzero(De,sizeof(De));
445: PetscMemzero(Ce,sizeof(Ce));
447: /* form element stiffness matrix */
448: BForm_DivT(Ae,el_coords,prop_eta);
449: BForm_Grad(Ge,el_coords);
450: BForm_ScaledMassMatrix(Ce,el_coords,prop_eta);
452: /* insert element matrix into global matrix */
453: DMDAGetElementEqnums_up(element,u_eqn,p_eqn);
454: MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
455: MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
456: MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
457: MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
458: }
459: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
460: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
462: DMSwarmRestoreField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
463: VecRestoreArrayRead(coords,&_coords);
465: return 0;
466: }
468: static PetscErrorCode AssembleStokes_RHS(Vec F,DM stokes_da,DM quadrature)
469: {
470: DM cda;
471: Vec coords;
472: const PetscScalar *_coords;
473: PetscInt u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
474: PetscInt p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
475: PetscInt nel,npe,eidx,i;
476: const PetscInt *element_list;
477: PetscScalar Fe[NODES_PER_EL*U_DOFS];
478: PetscScalar He[NODES_PER_EL*P_DOFS];
479: PetscScalar el_coords[NODES_PER_EL*NSD];
480: PetscScalar *q_rhs,*prop_fy;
481: Vec local_F;
482: PetscScalar *LA_F;
485: VecZeroEntries(F);
486: /* setup for coords */
487: DMGetCoordinateDM(stokes_da,&cda);
488: DMGetCoordinatesLocal(stokes_da,&coords);
489: VecGetArrayRead(coords,&_coords);
491: /* setup for coefficients */
492: DMSwarmGetField(quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
494: /* get access to the vector */
495: DMGetLocalVector(stokes_da,&local_F);
496: VecZeroEntries(local_F);
497: VecGetArray(local_F,&LA_F);
499: DMDAGetElements(stokes_da,&nel,&npe,&element_list);
500: for (eidx = 0; eidx < nel; eidx++) {
501: const PetscInt *element = &element_list[npe*eidx];
503: /* get coords for the element */
504: GetElementCoords(_coords,element,el_coords);
506: /* get coefficients for the element */
507: prop_fy = &q_rhs[GAUSS_POINTS * eidx];
509: /* initialise element stiffness matrix */
510: PetscMemzero(Fe,sizeof(Fe));
511: PetscMemzero(He,sizeof(He));
513: /* form element stiffness matrix */
514: LForm_MomentumRHS(Fe,el_coords,NULL,prop_fy);
516: /* insert element matrix into global matrix */
517: DMDAGetElementEqnums_up(element,u_eqn,p_eqn);
519: for (i=0; i<NODES_PER_EL*U_DOFS; i++) {
520: LA_F[ u_eqn[i] ] += Fe[i];
521: }
522: }
523: DMSwarmRestoreField(quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
524: VecRestoreArrayRead(coords,&_coords);
526: VecRestoreArray(local_F,&LA_F);
527: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
528: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
529: DMRestoreLocalVector(stokes_da,&local_F);
531: return 0;
532: }
534: PetscErrorCode DMSwarmPICInsertPointsCellwise(DM dm,DM dmc,PetscInt e,PetscInt npoints,PetscReal xi[],PetscBool proximity_initialization)
535: {
536: PetscInt dim,nel,npe,q,k,d,ncurr;
537: const PetscInt *element_list;
538: Vec coor;
539: const PetscScalar *_coor;
540: PetscReal **basis,*elcoor,*xp;
541: PetscReal *swarm_coor;
542: PetscInt *swarm_cellid;
545: DMGetDimension(dm,&dim);
546: DMDAGetElements(dmc,&nel,&npe,&element_list);
548: PetscMalloc1(dim*npoints,&xp);
549: PetscMalloc1(dim*npe,&elcoor);
550: PetscMalloc1(npoints,&basis);
551: for (q=0; q<npoints; q++) {
552: PetscMalloc1(npe,&basis[q]);
554: switch (dim) {
555: case 1:
556: basis[q][0] = 0.5*(1.0 - xi[dim*q+0]);
557: basis[q][1] = 0.5*(1.0 + xi[dim*q+0]);
558: break;
559: case 2:
560: basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]);
561: basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]);
562: basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]);
563: basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]);
564: break;
566: case 3:
567: basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
568: basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
569: basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
570: basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
571: basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
572: basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
573: basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
574: basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
575: break;
576: }
577: }
579: DMGetCoordinatesLocal(dmc,&coor);
580: VecGetArrayRead(coor,&_coor);
581: /* compute and store the coordinates for the new points */
582: {
583: const PetscInt *element = &element_list[npe*e];
585: for (k=0; k<npe; k++) {
586: for (d=0; d<dim; d++) {
587: elcoor[dim*k+d] = PetscRealPart(_coor[ dim*element[k] + d ]);
588: }
589: }
590: for (q=0; q<npoints; q++) {
591: for (d=0; d<dim; d++) {
592: xp[dim*q+d] = 0.0;
593: }
594: for (k=0; k<npe; k++) {
595: for (d=0; d<dim; d++) {
596: xp[dim*q+d] += basis[q][k] * elcoor[dim*k+d];
597: }
598: }
599: }
600: }
601: VecRestoreArrayRead(coor,&_coor);
602: DMDARestoreElements(dmc,&nel,&npe,&element_list);
604: DMSwarmGetLocalSize(dm,&ncurr);
605: DMSwarmAddNPoints(dm,npoints);
607: if (proximity_initialization) {
608: PetscInt *nnlist;
609: PetscReal *coor_q,*coor_qn;
610: PetscInt npoints_e,*plist_e;
612: DMSwarmSortGetPointsPerCell(dm,e,&npoints_e,&plist_e);
614: PetscMalloc1(npoints,&nnlist);
615: /* find nearest neighour points in this cell */
616: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
617: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
618: for (q=0; q<npoints; q++) {
619: PetscInt qn,nearest_neighour = -1;
620: PetscReal sep,min_sep = PETSC_MAX_REAL;
622: coor_q = &xp[dim*q];
623: for (qn=0; qn<npoints_e; qn++) {
624: coor_qn = &swarm_coor[dim*plist_e[qn]];
625: sep = 0.0;
626: for (d=0; d<dim; d++) {
627: sep += (coor_q[d]-coor_qn[d])*(coor_q[d]-coor_qn[d]);
628: }
629: if (sep < min_sep) {
630: nearest_neighour = plist_e[qn];
631: min_sep = sep;
632: }
633: }
635: nnlist[q] = nearest_neighour;
636: }
637: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
638: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
640: /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
641: for (q=0; q<npoints; q++) {
642: DMSwarmCopyPoint(dm,nnlist[q],ncurr+q);
643: }
644: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
645: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
646: for (q=0; q<npoints; q++) {
647: /* set the coordinates */
648: for (d=0; d<dim; d++) {
649: swarm_coor[dim*(ncurr+q)+d] = xp[dim*q+d];
650: }
651: /* set the cell index */
652: swarm_cellid[ncurr+q] = e;
653: }
654: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
655: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
657: PetscFree(plist_e);
658: PetscFree(nnlist);
659: } else {
660: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
661: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
662: for (q=0; q<npoints; q++) {
663: /* set the coordinates */
664: for (d=0; d<dim; d++) {
665: swarm_coor[dim*(ncurr+q)+d] = xp[dim*q+d];
666: }
667: /* set the cell index */
668: swarm_cellid[ncurr+q] = e;
669: }
670: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
671: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
672: }
674: PetscFree(xp);
675: PetscFree(elcoor);
676: for (q=0; q<npoints; q++) {
677: PetscFree(basis[q]);
678: }
679: PetscFree(basis);
680: return 0;
681: }
683: PetscErrorCode MaterialPoint_PopulateCell(DM dm_vp,DM dm_mpoint)
684: {
685: PetscInt _npe,_nel,e,nel;
686: const PetscInt *element;
687: DM dmc;
688: PetscQuadrature quadrature;
689: const PetscReal *xi;
690: PetscInt npoints_q,cnt,cnt_g;
693: DMDAGetElements(dm_vp,&_nel,&_npe,&element);
694: nel = _nel;
695: DMDARestoreElements(dm_vp,&_nel,&_npe,&element);
697: PetscDTGaussTensorQuadrature(2,1,4,-1.0,1.0,&quadrature);
698: PetscQuadratureGetData(quadrature,NULL,NULL,&npoints_q,&xi,NULL);
699: DMSwarmGetCellDM(dm_mpoint,&dmc);
701: DMSwarmSortGetAccess(dm_mpoint);
703: cnt = 0;
704: for (e=0; e<nel; e++) {
705: PetscInt npoints_per_cell;
707: DMSwarmSortGetNumberOfPointsPerCell(dm_mpoint,e,&npoints_per_cell);
709: if (npoints_per_cell < 12) {
710: DMSwarmPICInsertPointsCellwise(dm_mpoint,dm_vp,e,npoints_q,(PetscReal*)xi,PETSC_TRUE);
711: cnt++;
712: }
713: }
714: MPI_Allreduce(&cnt,&cnt_g,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
715: if (cnt_g > 0) {
716: PetscPrintf(PETSC_COMM_WORLD,".... ....pop cont: adjusted %D cells\n",cnt_g);
717: }
719: DMSwarmSortRestoreAccess(dm_mpoint);
720: PetscQuadratureDestroy(&quadrature);
721: return 0;
722: }
724: PetscErrorCode MaterialPoint_AdvectRK1(DM dm_vp,Vec vp,PetscReal dt,DM dm_mpoint)
725: {
726: Vec vp_l,coor_l;
727: const PetscScalar *LA_vp;
728: PetscInt i,p,e,npoints,nel,npe;
729: PetscInt *mpfield_cell;
730: PetscReal *mpfield_coor;
731: const PetscInt *element_list;
732: const PetscInt *element;
733: PetscScalar xi_p[NSD],Ni[NODES_PER_EL];
734: const PetscScalar *LA_coor;
735: PetscScalar dx[NSD];
738: DMGetCoordinatesLocal(dm_vp,&coor_l);
739: VecGetArrayRead(coor_l,&LA_coor);
741: DMGetLocalVector(dm_vp,&vp_l);
742: DMGlobalToLocalBegin(dm_vp,vp,INSERT_VALUES,vp_l);
743: DMGlobalToLocalEnd(dm_vp,vp,INSERT_VALUES,vp_l);
744: VecGetArrayRead(vp_l,&LA_vp);
746: DMDAGetElements(dm_vp,&nel,&npe,&element_list);
747: DMSwarmGetLocalSize(dm_mpoint,&npoints);
748: DMSwarmGetField(dm_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
749: DMSwarmGetField(dm_mpoint,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
750: for (p=0; p<npoints; p++) {
751: PetscReal *coor_p;
752: PetscScalar vel_n[NSD*NODES_PER_EL],vel_p[NSD];
753: const PetscScalar *x0;
754: const PetscScalar *x2;
756: e = mpfield_cell[p];
757: coor_p = &mpfield_coor[NSD*p];
758: element = &element_list[NODES_PER_EL*e];
760: /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
761: x0 = &LA_coor[NSD*element[0]];
762: x2 = &LA_coor[NSD*element[2]];
764: dx[0] = x2[0] - x0[0];
765: dx[1] = x2[1] - x0[1];
767: xi_p[0] = 2.0 * (coor_p[0] - x0[0])/dx[0] - 1.0;
768: xi_p[1] = 2.0 * (coor_p[1] - x0[1])/dx[1] - 1.0;
774: /* evaluate basis functions */
775: EvaluateBasis_Q1(xi_p,Ni);
777: /* get cell nodal velocities */
778: for (i=0; i<NODES_PER_EL; i++) {
779: PetscInt nid;
781: nid = element[i];
782: vel_n[NSD*i+0] = LA_vp[(NSD+1)*nid+0];
783: vel_n[NSD*i+1] = LA_vp[(NSD+1)*nid+1];
784: }
786: /* interpolate velocity */
787: vel_p[0] = vel_p[1] = 0.0;
788: for (i=0; i<NODES_PER_EL; i++) {
789: vel_p[0] += Ni[i] * vel_n[NSD*i+0];
790: vel_p[1] += Ni[i] * vel_n[NSD*i+1];
791: }
793: coor_p[0] += dt * PetscRealPart(vel_p[0]);
794: coor_p[1] += dt * PetscRealPart(vel_p[1]);
795: }
797: DMSwarmRestoreField(dm_mpoint,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
798: DMSwarmRestoreField(dm_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
799: DMDARestoreElements(dm_vp,&nel,&npe,&element_list);
800: VecRestoreArrayRead(vp_l,&LA_vp);
801: DMRestoreLocalVector(dm_vp,&vp_l);
802: VecRestoreArrayRead(coor_l,&LA_coor);
803: return 0;
804: }
806: PetscErrorCode MaterialPoint_Interpolate(DM dm,Vec eta_v,Vec rho_v,DM dm_quadrature)
807: {
808: Vec eta_l,rho_l;
809: PetscScalar *_eta_l,*_rho_l;
810: PetscInt nqp,npe,nel;
811: PetscScalar qp_xi[GAUSS_POINTS][NSD];
812: PetscScalar qp_weight[GAUSS_POINTS];
813: PetscInt q,k,e;
814: PetscScalar Ni[GAUSS_POINTS][NODES_PER_EL];
815: const PetscInt *element_list;
816: PetscReal *q_eta,*q_rhs;
819: /* define quadrature rule */
820: CreateGaussQuadrature(&nqp,qp_xi,qp_weight);
821: for (q=0; q<nqp; q++) {
822: EvaluateBasis_Q1(qp_xi[q],Ni[q]);
823: }
825: DMGetLocalVector(dm,&eta_l);
826: DMGetLocalVector(dm,&rho_l);
828: DMGlobalToLocalBegin(dm,eta_v,INSERT_VALUES,eta_l);
829: DMGlobalToLocalEnd(dm,eta_v,INSERT_VALUES,eta_l);
830: DMGlobalToLocalBegin(dm,rho_v,INSERT_VALUES,rho_l);
831: DMGlobalToLocalEnd(dm,rho_v,INSERT_VALUES,rho_l);
833: VecGetArray(eta_l,&_eta_l);
834: VecGetArray(rho_l,&_rho_l);
836: DMSwarmGetField(dm_quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
837: DMSwarmGetField(dm_quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
839: DMDAGetElements(dm,&nel,&npe,&element_list);
840: for (e=0; e<nel; e++) {
841: PetscScalar eta_field_e[NODES_PER_EL];
842: PetscScalar rho_field_e[NODES_PER_EL];
843: const PetscInt *element = &element_list[4*e];
845: for (k=0; k<NODES_PER_EL; k++) {
846: eta_field_e[k] = _eta_l[ element[k] ];
847: rho_field_e[k] = _rho_l[ element[k] ];
848: }
850: for (q=0; q<nqp; q++) {
851: PetscScalar eta_q,rho_q;
853: eta_q = rho_q = 0.0;
854: for (k=0; k<NODES_PER_EL; k++) {
855: eta_q += Ni[q][k] * eta_field_e[k];
856: rho_q += Ni[q][k] * rho_field_e[k];
857: }
859: q_eta[nqp*e+q] = PetscRealPart(eta_q);
860: q_rhs[nqp*e+q] = PetscRealPart(rho_q);
861: }
862: }
863: DMDARestoreElements(dm,&nel,&npe,&element_list);
865: DMSwarmRestoreField(dm_quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
866: DMSwarmRestoreField(dm_quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
868: VecRestoreArray(rho_l,&_rho_l);
869: VecRestoreArray(eta_l,&_eta_l);
870: DMRestoreLocalVector(dm,&rho_l);
871: DMRestoreLocalVector(dm,&eta_l);
872: return 0;
873: }
875: static PetscErrorCode SolveTimeDepStokes(PetscInt mx,PetscInt my)
876: {
877: DM dm_stokes,dm_coeff;
878: PetscInt u_dof,p_dof,dof,stencil_width;
879: Mat A,B;
880: PetscInt nel_local;
881: Vec eta_v,rho_v;
882: Vec f,X;
883: KSP ksp;
884: PC pc;
885: char filename[PETSC_MAX_PATH_LEN];
886: DM dms_quadrature,dms_mpoint;
887: PetscInt nel,npe,npoints;
888: const PetscInt *element_list;
889: PetscInt tk,nt,dump_freq;
890: PetscReal dt,dt_max = 0.0;
891: PetscReal vx[2],vy[2],max_v = 0.0,max_v_step,dh;
892: const char *fieldnames[] = { "eta" , "rho" };
893: Vec *pfields;
894: PetscInt ppcell = 1;
895: PetscReal time,delta_eta = 1.0;
896: PetscBool randomize_coords = PETSC_FALSE;
897: PetscReal randomize_fac = 0.25;
898: PetscBool no_view = PETSC_FALSE;
899: PetscBool isbddc;
902: /*
903: Generate the DMDA for the velocity and pressure spaces.
904: We use Q1 elements for both fields.
905: The Q1 FE basis on a regular mesh has a 9-point stencil (DMDA_STENCIL_BOX)
906: The number of nodes in each direction is mx+1, my+1
907: */
908: u_dof = U_DOFS; /* Vx, Vy - velocities */
909: p_dof = P_DOFS; /* p - pressure */
910: dof = u_dof + p_dof;
911: stencil_width = 1;
912: 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);
913: DMDASetElementType(dm_stokes,DMDA_ELEMENT_Q1);
914: DMSetMatType(dm_stokes,MATAIJ);
915: DMSetFromOptions(dm_stokes);
916: DMSetUp(dm_stokes);
917: DMDASetFieldName(dm_stokes,0,"ux");
918: DMDASetFieldName(dm_stokes,1,"uy");
919: DMDASetFieldName(dm_stokes,2,"p");
921: /* unit box [0,0.9142] x [0,1] */
922: DMDASetUniformCoordinates(dm_stokes,0.0,0.9142,0.0,1.0,0.,0.);
923: dh = 1.0/((PetscReal)(mx));
925: /* Get local number of elements */
926: {
927: DMDAGetElements(dm_stokes,&nel,&npe,&element_list);
929: nel_local = nel;
931: DMDARestoreElements(dm_stokes,&nel,&npe,&element_list);
932: }
934: /* Create DMDA for representing scalar fields */
935: DMDACreateCompatibleDMDA(dm_stokes,1,&dm_coeff);
937: /* Create the swarm for storing quadrature point values */
938: DMCreate(PETSC_COMM_WORLD,&dms_quadrature);
939: DMSetType(dms_quadrature,DMSWARM);
940: DMSetDimension(dms_quadrature,2);
942: /* Register fields for viscosity and density on the quadrature points */
943: DMSwarmRegisterPetscDatatypeField(dms_quadrature,"eta_q",1,PETSC_REAL);
944: DMSwarmRegisterPetscDatatypeField(dms_quadrature,"rho_q",1,PETSC_REAL);
945: DMSwarmFinalizeFieldRegister(dms_quadrature);
946: DMSwarmSetLocalSizes(dms_quadrature,nel_local * GAUSS_POINTS,0);
948: /* Create the material point swarm */
949: DMCreate(PETSC_COMM_WORLD,&dms_mpoint);
950: DMSetType(dms_mpoint,DMSWARM);
951: DMSetDimension(dms_mpoint,2);
953: /* Configure the material point swarm to be of type Particle-In-Cell */
954: DMSwarmSetType(dms_mpoint,DMSWARM_PIC);
956: /*
957: Specify the DM to use for point location and projections
958: within the context of a PIC scheme
959: */
960: DMSwarmSetCellDM(dms_mpoint,dm_coeff);
962: /* Register fields for viscosity and density */
963: DMSwarmRegisterPetscDatatypeField(dms_mpoint,"eta",1,PETSC_REAL);
964: DMSwarmRegisterPetscDatatypeField(dms_mpoint,"rho",1,PETSC_REAL);
965: DMSwarmFinalizeFieldRegister(dms_mpoint);
967: PetscOptionsGetInt(NULL,NULL,"-ppcell",&ppcell,NULL);
968: DMSwarmSetLocalSizes(dms_mpoint,nel_local * ppcell,100);
970: /*
971: Layout the material points in space using the cell DM.
972: Particle coordinates are defined by cell wise using different methods.
973: - DMSWARMPIC_LAYOUT_GAUSS defines particles coordinates at the positions
974: corresponding to a Gauss quadrature rule with
975: ppcell points in each direction.
976: - DMSWARMPIC_LAYOUT_REGULAR defines particle coordinates at the centoid of
977: ppcell x ppcell quadralaterals defined within the
978: reference element.
979: - DMSWARMPIC_LAYOUT_SUBDIVISION defines particles coordinates at the centroid
980: of each quadralateral obtained by sub-dividing
981: the reference element cell ppcell times.
982: */
983: DMSwarmInsertPointsUsingCellDM(dms_mpoint,DMSWARMPIC_LAYOUT_SUBDIVISION,ppcell);
985: /*
986: Defne a high resolution layer of material points across the material interface
987: */
988: {
989: PetscInt npoints_dir_x[2];
990: PetscReal min[2],max[2];
992: npoints_dir_x[0] = (PetscInt)(0.9142/(0.05*dh));
993: npoints_dir_x[1] = (PetscInt)((0.25-0.15)/(0.05*dh));
994: min[0] = 0.0; max[0] = 0.9142;
995: min[1] = 0.05; max[1] = 0.35;
996: DMSwarmSetPointsUniformCoordinates(dms_mpoint,min,max,npoints_dir_x,ADD_VALUES);
997: }
999: /*
1000: Define a high resolution layer of material points near the surface of the domain
1001: to deal with weakly compressible Q1-Q1 elements. These elements "self compact"
1002: when applied to buouyancy driven flow. The error in div(u) is O(h).
1003: */
1004: {
1005: PetscInt npoints_dir_x[2];
1006: PetscReal min[2],max[2];
1008: npoints_dir_x[0] = (PetscInt)(0.9142/(0.25*dh));
1009: npoints_dir_x[1] = (PetscInt)(3.0*dh/(0.25*dh));
1010: min[0] = 0.0; max[0] = 0.9142;
1011: min[1] = 1.0 - 3.0*dh; max[1] = 1.0-0.0001;
1012: DMSwarmSetPointsUniformCoordinates(dms_mpoint,min,max,npoints_dir_x,ADD_VALUES);
1013: }
1015: DMView(dms_mpoint,PETSC_VIEWER_STDOUT_WORLD);
1017: /* Define initial material properties on each particle in the material point swarm */
1018: PetscOptionsGetReal(NULL,NULL,"-delta_eta",&delta_eta,NULL);
1019: PetscOptionsGetBool(NULL,NULL,"-randomize_coords",&randomize_coords,NULL);
1020: PetscOptionsGetReal(NULL,NULL,"-randomize_fac",&randomize_fac,NULL);
1022: {
1023: PetscReal *array_x,*array_e,*array_r;
1024: PetscInt p;
1025: PetscRandom r;
1026: PetscMPIInt rank;
1028: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1030: PetscRandomCreate(PETSC_COMM_SELF,&r);
1031: PetscRandomSetInterval(r,-randomize_fac*dh,randomize_fac*dh);
1032: PetscRandomSetSeed(r,(unsigned long)rank);
1033: PetscRandomSeed(r);
1035: DMDAGetElements(dm_stokes,&nel,&npe,&element_list);
1037: /*
1038: Fetch the registered data from the material point DMSwarm.
1039: The fields "eta" and "rho" were registered by this example.
1040: The field identified by the the variable DMSwarmPICField_coor
1041: was registered by the DMSwarm implementation when the function
1042: DMSwarmSetType(dms_mpoint,DMSWARM_PIC)
1043: was called. The returned array defines the coordinates of each
1044: material point in the point swarm.
1045: */
1046: DMSwarmGetField(dms_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&array_x);
1047: DMSwarmGetField(dms_mpoint,"eta", NULL,NULL,(void**)&array_e);
1048: DMSwarmGetField(dms_mpoint,"rho", NULL,NULL,(void**)&array_r);
1050: DMSwarmGetLocalSize(dms_mpoint,&npoints);
1051: for (p = 0; p < npoints; p++) {
1052: PetscReal x_p[2],rr[2];
1054: if (randomize_coords) {
1055: PetscRandomGetValueReal(r,&rr[0]);
1056: PetscRandomGetValueReal(r,&rr[1]);
1057: array_x[2*p + 0] += rr[0];
1058: array_x[2*p + 1] += rr[1];
1059: }
1061: /* Get the coordinates of point, p */
1062: x_p[0] = array_x[2*p + 0];
1063: x_p[1] = array_x[2*p + 1];
1065: if (x_p[1] < (0.2 + 0.02*PetscCosReal(PETSC_PI*x_p[0]/0.9142))) {
1066: /* Material properties below the interface */
1067: array_e[p] = 1.0 * (1.0/delta_eta);
1068: array_r[p] = 0.0;
1069: } else {
1070: /* Material properties above the interface */
1071: array_e[p] = 1.0;
1072: array_r[p] = 1.0;
1073: }
1074: }
1076: /*
1077: Restore the fetched data fields from the material point DMSwarm.
1078: Calling the Restore function invalidates the points array_r, array_e, array_x
1079: by setting them to NULL.
1080: */
1081: DMSwarmRestoreField(dms_mpoint,"rho",NULL,NULL,(void**)&array_r);
1082: DMSwarmRestoreField(dms_mpoint,"eta",NULL,NULL,(void**)&array_e);
1083: DMSwarmRestoreField(dms_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&array_x);
1085: DMDARestoreElements(dm_stokes,&nel,&npe,&element_list);
1086: PetscRandomDestroy(&r);
1087: }
1089: /*
1090: If the particle coordinates where randomly shifted, they may have crossed into another
1091: element, or into another sub-domain. To account for this we call the Migrate function.
1092: */
1093: if (randomize_coords) {
1094: DMSwarmMigrate(dms_mpoint,PETSC_TRUE);
1095: }
1097: PetscOptionsGetBool(NULL,NULL,"-no_view",&no_view,NULL);
1098: if (!no_view) {
1099: DMSwarmViewXDMF(dms_mpoint,"ic_coeff_dms.xmf");
1100: }
1102: /* project the swarm properties */
1103: DMSwarmProjectFields(dms_mpoint,2,fieldnames,&pfields,PETSC_FALSE);
1104: eta_v = pfields[0];
1105: rho_v = pfields[1];
1106: PetscObjectSetName((PetscObject)eta_v,"eta");
1107: PetscObjectSetName((PetscObject)rho_v,"rho");
1108: MaterialPoint_Interpolate(dm_coeff,eta_v,rho_v,dms_quadrature);
1110: /* view projected coefficients eta and rho */
1111: if (!no_view) {
1112: PetscViewer viewer;
1114: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1115: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1116: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1117: PetscViewerFileSetName(viewer,"ic_coeff_dmda.vts");
1118: VecView(eta_v,viewer);
1119: VecView(rho_v,viewer);
1120: PetscViewerDestroy(&viewer);
1121: }
1123: DMCreateMatrix(dm_stokes,&A);
1124: DMCreateMatrix(dm_stokes,&B);
1125: DMCreateGlobalVector(dm_stokes,&f);
1126: DMCreateGlobalVector(dm_stokes,&X);
1128: AssembleStokes_A(A,dm_stokes,dms_quadrature);
1129: AssembleStokes_PC(B,dm_stokes,dms_quadrature);
1130: AssembleStokes_RHS(f,dm_stokes,dms_quadrature);
1132: DMDAApplyBoundaryConditions(dm_stokes,A,f);
1133: DMDAApplyBoundaryConditions(dm_stokes,B,NULL);
1135: KSPCreate(PETSC_COMM_WORLD,&ksp);
1136: KSPSetOptionsPrefix(ksp,"stokes_");
1137: KSPSetDM(ksp,dm_stokes);
1138: KSPSetDMActive(ksp,PETSC_FALSE);
1139: KSPSetOperators(ksp,A,B);
1140: KSPSetFromOptions(ksp);
1141: KSPGetPC(ksp,&pc);
1142: PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);
1143: if (isbddc) {
1144: KSPSetOperators(ksp,A,A);
1145: }
1147: /* Define u-v-p indices for fieldsplit */
1148: {
1149: PC pc;
1150: const PetscInt ufields[] = {0,1},pfields[1] = {2};
1152: KSPGetPC(ksp,&pc);
1153: PCFieldSplitSetBlockSize(pc,3);
1154: PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1155: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1156: }
1158: /* If using a fieldsplit preconditioner, attach a DMDA to the velocity split so that geometric multigrid can be used */
1159: {
1160: PC pc,pc_u;
1161: KSP *sub_ksp,ksp_u;
1162: PetscInt nsplits;
1163: DM dm_u;
1164: PetscBool is_pcfs;
1166: KSPGetPC(ksp,&pc);
1168: is_pcfs = PETSC_FALSE;
1169: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&is_pcfs);
1171: if (is_pcfs) {
1172: KSPSetUp(ksp);
1173: KSPGetPC(ksp,&pc);
1174: PCFieldSplitGetSubKSP(pc,&nsplits,&sub_ksp);
1175: ksp_u = sub_ksp[0];
1176: PetscFree(sub_ksp);
1178: if (nsplits == 2) {
1179: DMDACreateCompatibleDMDA(dm_stokes,2,&dm_u);
1181: KSPSetDM(ksp_u,dm_u);
1182: KSPSetDMActive(ksp_u,PETSC_FALSE);
1183: DMDestroy(&dm_u);
1185: /* enforce galerkin coarse grids be used */
1186: KSPGetPC(ksp_u,&pc_u);
1187: PCMGSetGalerkin(pc_u,PC_MG_GALERKIN_PMAT);
1188: }
1189: }
1190: }
1192: dump_freq = 10;
1193: PetscOptionsGetInt(NULL,NULL,"-dump_freq",&dump_freq,NULL);
1194: nt = 10;
1195: PetscOptionsGetInt(NULL,NULL,"-nt",&nt,NULL);
1196: time = 0.0;
1197: for (tk=1; tk<=nt; tk++) {
1199: PetscPrintf(PETSC_COMM_WORLD,".... assemble\n");
1200: AssembleStokes_A(A,dm_stokes,dms_quadrature);
1201: AssembleStokes_PC(B,dm_stokes,dms_quadrature);
1202: AssembleStokes_RHS(f,dm_stokes,dms_quadrature);
1204: PetscPrintf(PETSC_COMM_WORLD,".... bc imposition\n");
1205: DMDAApplyBoundaryConditions(dm_stokes,A,f);
1206: DMDAApplyBoundaryConditions(dm_stokes,B,NULL);
1208: PetscPrintf(PETSC_COMM_WORLD,".... solve\n");
1209: KSPSetOperators(ksp,A, isbddc ? A : B);
1210: KSPSolve(ksp,f,X);
1212: VecStrideMax(X,0,NULL,&vx[1]);
1213: VecStrideMax(X,1,NULL,&vy[1]);
1214: VecStrideMin(X,0,NULL,&vx[0]);
1215: VecStrideMin(X,1,NULL,&vy[0]);
1217: max_v_step = PetscMax(vx[0],vx[1]);
1218: max_v_step = PetscMax(max_v_step,vy[0]);
1219: max_v_step = PetscMax(max_v_step,vy[1]);
1220: max_v = PetscMax(max_v,max_v_step);
1222: dt_max = 2.0;
1223: dt = 0.5 * (dh / max_v_step);
1224: 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);
1225: dt = PetscMin(dt_max,dt);
1227: /* advect */
1228: PetscPrintf(PETSC_COMM_WORLD,".... advect\n");
1229: MaterialPoint_AdvectRK1(dm_stokes,X,dt,dms_mpoint);
1231: /* migrate */
1232: PetscPrintf(PETSC_COMM_WORLD,".... migrate\n");
1233: DMSwarmMigrate(dms_mpoint,PETSC_TRUE);
1235: /* update cell population */
1236: PetscPrintf(PETSC_COMM_WORLD,".... populate cells\n");
1237: MaterialPoint_PopulateCell(dm_stokes,dms_mpoint);
1239: /* update coefficients on quadrature points */
1240: PetscPrintf(PETSC_COMM_WORLD,".... project\n");
1241: DMSwarmProjectFields(dms_mpoint,2,fieldnames,&pfields,PETSC_TRUE);
1242: eta_v = pfields[0];
1243: rho_v = pfields[1];
1244: PetscPrintf(PETSC_COMM_WORLD,".... interp\n");
1245: MaterialPoint_Interpolate(dm_coeff,eta_v,rho_v,dms_quadrature);
1247: if (tk%dump_freq == 0) {
1248: PetscViewer viewer;
1250: PetscPrintf(PETSC_COMM_WORLD,".... write XDMF, VTS\n");
1251: PetscSNPrintf(filename,PETSC_MAX_PATH_LEN-1,"step%.4D_coeff_dms.xmf",tk);
1252: DMSwarmViewXDMF(dms_mpoint,filename);
1254: PetscSNPrintf(filename,PETSC_MAX_PATH_LEN-1,"step%.4D_vp_dm.vts",tk);
1255: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1256: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1257: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1258: PetscViewerFileSetName(viewer,filename);
1259: VecView(X,viewer);
1260: PetscViewerDestroy(&viewer);
1261: }
1262: time += dt;
1263: PetscPrintf(PETSC_COMM_WORLD,"step %D : time %1.2e \n",tk,(double)time);
1264: }
1266: KSPDestroy(&ksp);
1267: VecDestroy(&X);
1268: VecDestroy(&f);
1269: MatDestroy(&A);
1270: MatDestroy(&B);
1271: VecDestroy(&eta_v);
1272: VecDestroy(&rho_v);
1273: PetscFree(pfields);
1275: DMDestroy(&dms_mpoint);
1276: DMDestroy(&dms_quadrature);
1277: DMDestroy(&dm_coeff);
1278: DMDestroy(&dm_stokes);
1279: return 0;
1280: }
1282: /*
1283: <sequential run>
1284: ./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
1285: */
1286: int main(int argc,char **args)
1287: {
1288: PetscInt mx,my;
1289: PetscBool set = PETSC_FALSE;
1291: PetscInitialize(&argc,&args,(char*)0,help);
1292: mx = my = 10;
1293: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1294: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1295: PetscOptionsGetInt(NULL,NULL,"-mxy",&mx,&set);
1296: if (set) {
1297: my = mx;
1298: }
1299: SolveTimeDepStokes(mx,my);
1300: PetscFinalize();
1301: return 0;
1302: }
1304: /* -------------------------- helpers for boundary conditions -------------------------------- */
1305: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
1306: {
1307: DM cda;
1308: Vec coords;
1309: PetscInt si,sj,nx,ny,i,j;
1310: PetscInt M,N;
1311: DMDACoor2d **_coords;
1312: const PetscInt *g_idx;
1313: PetscInt *bc_global_ids;
1314: PetscScalar *bc_vals;
1315: PetscInt nbcs;
1316: PetscInt n_dofs;
1317: ISLocalToGlobalMapping ltogm;
1320: DMGetLocalToGlobalMapping(da,<ogm);
1321: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1323: DMGetCoordinateDM(da,&cda);
1324: DMGetCoordinatesLocal(da,&coords);
1325: DMDAVecGetArray(cda,coords,&_coords);
1326: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1327: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1329: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1330: PetscMalloc1(ny*n_dofs,&bc_vals);
1332: /* init the entries to -1 so VecSetValues will ignore them */
1333: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1335: i = nx-1;
1336: for (j = 0; j < ny; j++) {
1337: PetscInt local_id;
1339: local_id = i+j*nx;
1341: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1343: bc_vals[j] = 0.0;
1344: }
1345: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1346: nbcs = 0;
1347: if ((si+nx) == (M)) nbcs = ny;
1349: if (b) {
1350: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1351: VecAssemblyBegin(b);
1352: VecAssemblyEnd(b);
1353: }
1354: if (A) {
1355: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1356: }
1358: PetscFree(bc_vals);
1359: PetscFree(bc_global_ids);
1361: DMDAVecRestoreArray(cda,coords,&_coords);
1362: return 0;
1363: }
1365: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
1366: {
1367: DM cda;
1368: Vec coords;
1369: PetscInt si,sj,nx,ny,i,j;
1370: PetscInt M,N;
1371: DMDACoor2d **_coords;
1372: const PetscInt *g_idx;
1373: PetscInt *bc_global_ids;
1374: PetscScalar *bc_vals;
1375: PetscInt nbcs;
1376: PetscInt n_dofs;
1377: ISLocalToGlobalMapping ltogm;
1380: DMGetLocalToGlobalMapping(da,<ogm);
1381: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1383: DMGetCoordinateDM(da,&cda);
1384: DMGetCoordinatesLocal(da,&coords);
1385: DMDAVecGetArray(cda,coords,&_coords);
1386: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1387: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1389: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1390: PetscMalloc1(ny*n_dofs,&bc_vals);
1392: /* init the entries to -1 so VecSetValues will ignore them */
1393: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1395: i = 0;
1396: for (j = 0; j < ny; j++) {
1397: PetscInt local_id;
1399: local_id = i+j*nx;
1401: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1403: bc_vals[j] = 0.0;
1404: }
1405: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1406: nbcs = 0;
1407: if (si == 0) nbcs = ny;
1409: if (b) {
1410: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1411: VecAssemblyBegin(b);
1412: VecAssemblyEnd(b);
1413: }
1415: if (A) {
1416: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1417: }
1419: PetscFree(bc_vals);
1420: PetscFree(bc_global_ids);
1422: DMDAVecRestoreArray(cda,coords,&_coords);
1423: return 0;
1424: }
1426: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
1427: {
1428: DM cda;
1429: Vec coords;
1430: PetscInt si,sj,nx,ny,i,j;
1431: PetscInt M,N;
1432: DMDACoor2d **_coords;
1433: const PetscInt *g_idx;
1434: PetscInt *bc_global_ids;
1435: PetscScalar *bc_vals;
1436: PetscInt nbcs;
1437: PetscInt n_dofs;
1438: ISLocalToGlobalMapping ltogm;
1441: DMGetLocalToGlobalMapping(da,<ogm);
1442: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1444: DMGetCoordinateDM(da,&cda);
1445: DMGetCoordinatesLocal(da,&coords);
1446: DMDAVecGetArray(cda,coords,&_coords);
1447: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1448: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1450: PetscMalloc1(nx,&bc_global_ids);
1451: PetscMalloc1(nx,&bc_vals);
1453: /* init the entries to -1 so VecSetValues will ignore them */
1454: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1456: j = ny-1;
1457: for (i = 0; i < nx; i++) {
1458: PetscInt local_id;
1460: local_id = i+j*nx;
1462: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1464: bc_vals[i] = 0.0;
1465: }
1466: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1467: nbcs = 0;
1468: if ((sj+ny) == (N)) nbcs = nx;
1470: if (b) {
1471: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1472: VecAssemblyBegin(b);
1473: VecAssemblyEnd(b);
1474: }
1475: if (A) {
1476: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,NULL,NULL);
1477: }
1479: PetscFree(bc_vals);
1480: PetscFree(bc_global_ids);
1482: DMDAVecRestoreArray(cda,coords,&_coords);
1483: return 0;
1484: }
1486: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
1487: {
1488: DM cda;
1489: Vec coords;
1490: PetscInt si,sj,nx,ny,i,j;
1491: PetscInt M,N;
1492: DMDACoor2d **_coords;
1493: const PetscInt *g_idx;
1494: PetscInt *bc_global_ids;
1495: PetscScalar *bc_vals;
1496: PetscInt nbcs;
1497: PetscInt n_dofs;
1498: ISLocalToGlobalMapping ltogm;
1501: DMGetLocalToGlobalMapping(da,<ogm);
1502: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1504: DMGetCoordinateDM(da,&cda);
1505: DMGetCoordinatesLocal(da,&coords);
1506: DMDAVecGetArray(cda,coords,&_coords);
1507: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1508: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1510: PetscMalloc1(nx,&bc_global_ids);
1511: PetscMalloc1(nx,&bc_vals);
1513: /* init the entries to -1 so VecSetValues will ignore them */
1514: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1516: j = 0;
1517: for (i = 0; i < nx; i++) {
1518: PetscInt local_id;
1520: local_id = i+j*nx;
1522: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1524: bc_vals[i] = 0.0;
1525: }
1526: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1527: nbcs = 0;
1528: if (sj == 0) nbcs = nx;
1530: if (b) {
1531: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1532: VecAssemblyBegin(b);
1533: VecAssemblyEnd(b);
1534: }
1535: if (A) {
1536: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1537: }
1539: PetscFree(bc_vals);
1540: PetscFree(bc_global_ids);
1542: DMDAVecRestoreArray(cda,coords,&_coords);
1543: return 0;
1544: }
1546: /*
1547: Impose free slip boundary conditions on the left/right faces: u_i n_i = 0, tau_{ij} t_j = 0
1548: Impose no slip boundray conditions on the top/bottom faces: u_i n_i = 0, u_i t_i = 0
1549: */
1550: static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes,Mat A,Vec f)
1551: {
1553: BCApplyZero_NORTH(dm_stokes,0,A,f);
1554: BCApplyZero_NORTH(dm_stokes,1,A,f);
1555: BCApplyZero_EAST(dm_stokes,0,A,f);
1556: BCApplyZero_SOUTH(dm_stokes,0,A,f);
1557: BCApplyZero_SOUTH(dm_stokes,1,A,f);
1558: BCApplyZero_WEST(dm_stokes,0,A,f);
1559: return 0;
1560: }
1562: /*TEST
1564: test:
1565: suffix: 1
1566: args: -no_view
1567: requires: !complex double
1568: filter: grep -v atomic
1569: filter_output: grep -v atomic
1570: test:
1571: suffix: 1_matis
1572: requires: !complex double
1573: args: -no_view -dm_mat_type is
1574: filter: grep -v atomic
1575: filter_output: grep -v atomic
1576: testset:
1577: nsize: 4
1578: requires: !complex double
1579: 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
1580: filter: grep -v atomic
1581: filter_output: grep -v atomic
1582: test:
1583: suffix: fetidp
1584: args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
1585: test:
1586: suffix: fetidp_lumped
1587: 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
1588: test:
1589: suffix: fetidp_saddlepoint
1590: 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
1591: test:
1592: suffix: fetidp_saddlepoint_lumped
1593: 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
1594: TEST*/