Actual source code: ex70.c
petsc-3.10.5 2019-03-28
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/examples/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
45: static void EvaluateBasis_Q1(PetscScalar _xi[],PetscScalar N[])
46: {
47: PetscScalar xi = _xi[0];
48: PetscScalar eta = _xi[1];
50: N[0] = 0.25*(1.0-xi)*(1.0-eta);
51: N[1] = 0.25*(1.0+xi)*(1.0-eta);
52: N[2] = 0.25*(1.0+xi)*(1.0+eta);
53: N[3] = 0.25*(1.0-xi)*(1.0+eta);
54: }
56: static void EvaluateBasisDerivatives_Q1(PetscScalar _xi[],PetscScalar dN[][NODES_PER_EL])
57: {
58: PetscScalar xi = _xi[0];
59: PetscScalar eta = _xi[1];
61: dN[0][0] = -0.25*(1.0-eta);
62: dN[0][1] = 0.25*(1.0-eta);
63: dN[0][2] = 0.25*(1.0+eta);
64: dN[0][3] = -0.25*(1.0+eta);
66: dN[1][0] = -0.25*(1.0-xi);
67: dN[1][1] = -0.25*(1.0+xi);
68: dN[1][2] = 0.25*(1.0+xi);
69: dN[1][3] = 0.25*(1.0-xi);
70: }
72: static void EvaluateDerivatives(PetscScalar dN[][NODES_PER_EL],PetscScalar dNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar *det_J)
73: {
74: PetscScalar J00,J01,J10,J11,J;
75: PetscScalar iJ00,iJ01,iJ10,iJ11;
76: PetscInt i;
78: J00 = J01 = J10 = J11 = 0.0;
79: for (i = 0; i < NODES_PER_EL; i++) {
80: PetscScalar cx = coords[2*i];
81: PetscScalar cy = coords[2*i+1];
83: J00 += dN[0][i]*cx; /* J_xx = dx/dxi */
84: J01 += dN[0][i]*cy; /* J_xy = dy/dxi */
85: J10 += dN[1][i]*cx; /* J_yx = dx/deta */
86: J11 += dN[1][i]*cy; /* J_yy = dy/deta */
87: }
88: J = (J00*J11)-(J01*J10);
90: iJ00 = J11/J;
91: iJ01 = -J01/J;
92: iJ10 = -J10/J;
93: iJ11 = J00/J;
95: for (i = 0; i < NODES_PER_EL; i++) {
96: dNx[0][i] = dN[0][i]*iJ00+dN[1][i]*iJ01;
97: dNx[1][i] = dN[0][i]*iJ10+dN[1][i]*iJ11;
98: }
100: if (det_J) *det_J = J;
101: }
103: static void CreateGaussQuadrature(PetscInt *ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])
104: {
105: *ngp = 4;
106: gp_xi[0][0] = -0.57735026919; gp_xi[0][1] = -0.57735026919;
107: gp_xi[1][0] = -0.57735026919; gp_xi[1][1] = 0.57735026919;
108: gp_xi[2][0] = 0.57735026919; gp_xi[2][1] = 0.57735026919;
109: gp_xi[3][0] = 0.57735026919; gp_xi[3][1] = -0.57735026919;
110: gp_weight[0] = 1.0;
111: gp_weight[1] = 1.0;
112: gp_weight[2] = 1.0;
113: gp_weight[3] = 1.0;
114: }
116: static PetscErrorCode DMDAGetElementEqnums_up(const PetscInt element[],PetscInt s_u[],PetscInt s_p[])
117: {
118: PetscInt i;
120: for (i=0; i<NODES_PER_EL; i++) {
121: /* velocity */
122: s_u[NSD*i+0] = 3*element[i];
123: s_u[NSD*i+1] = 3*element[i]+1;
124: /* pressure */
125: s_p[i] = 3*element[i]+2;
126: }
127: return(0);
128: }
130: 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)
131: {
132: PetscInt ij,r,c,nc;
134: nc = u_NPE*u_dof;
135: r = w_dof*wi+wd;
136: c = u_dof*ui+ud;
137: ij = r*nc+c;
138: return(ij);
139: }
141: static void BForm_DivT(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
142: {
143: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
144: PetscScalar GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
145: PetscScalar J_p,tildeD[3];
146: PetscScalar B[3][U_DOFS*NODES_PER_EL];
147: PetscInt p,i,j,k,ngp;
149: /* define quadrature rule */
150: CreateGaussQuadrature(&ngp,gp_xi,gp_weight);
152: /* evaluate bilinear form */
153: for (p = 0; p < ngp; p++) {
154: EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
155: EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
157: for (i = 0; i < NODES_PER_EL; i++) {
158: PetscScalar d_dx_i = GNx_p[0][i];
159: PetscScalar d_dy_i = GNx_p[1][i];
161: B[0][2*i] = d_dx_i;B[0][2*i+1] = 0.0;
162: B[1][2*i] = 0.0;B[1][2*i+1] = d_dy_i;
163: B[2][2*i] = d_dy_i;B[2][2*i+1] = d_dx_i;
164: }
166: tildeD[0] = 2.0*gp_weight[p]*J_p*eta[p];
167: tildeD[1] = 2.0*gp_weight[p]*J_p*eta[p];
168: tildeD[2] = gp_weight[p]*J_p*eta[p];
170: /* form Bt tildeD B */
171: /*
172: Ke_ij = Bt_ik . D_kl . B_lj
173: = B_ki . D_kl . B_lj
174: = B_ki . D_kk . B_kj
175: */
176: for (i = 0; i < 8; i++) {
177: for (j = 0; j < 8; j++) {
178: for (k = 0; k < 3; k++) { /* Note D is diagonal for stokes */
179: Ke[i+8*j] += B[k][i]*tildeD[k]*B[k][j];
180: }
181: }
182: }
183: }
184: }
186: static void BForm_Grad(PetscScalar Ke[],PetscScalar coords[])
187: {
188: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
189: PetscScalar Ni_p[NODES_PER_EL],GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
190: PetscScalar J_p,fac;
191: PetscInt p,i,j,di,ngp;
193: /* define quadrature rule */
194: CreateGaussQuadrature(&ngp,gp_xi,gp_weight);
196: /* evaluate bilinear form */
197: for (p = 0; p < ngp; p++) {
198: EvaluateBasis_Q1(gp_xi[p],Ni_p);
199: EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
200: EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
201: fac = gp_weight[p]*J_p;
203: for (i = 0; i < NODES_PER_EL; i++) { /* u nodes */
204: for (di = 0; di < NSD; di++) { /* u dofs */
205: for (j = 0; j < 4; j++) { /* p nodes, p dofs = 1 (ie no loop) */
206: PetscInt IJ;
207: IJ = map_wIwDI_uJuDJ(i,di,NODES_PER_EL,2,j,0,NODES_PER_EL,1);
209: Ke[IJ] -= GNx_p[di][i]*Ni_p[j]*fac;
210: }
211: }
212: }
213: }
214: }
216: static void BForm_Div(PetscScalar De[],PetscScalar coords[])
217: {
218: PetscScalar Ge[U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL];
219: PetscInt i,j,nr_g,nc_g;
221: PetscMemzero(Ge,sizeof(PetscScalar)*U_DOFS*NODES_PER_EL*P_DOFS*NODES_PER_EL);
222: BForm_Grad(Ge,coords);
224: nr_g = U_DOFS*NODES_PER_EL;
225: nc_g = P_DOFS*NODES_PER_EL;
227: for (i = 0; i < nr_g; i++) {
228: for (j = 0; j < nc_g; j++) {
229: De[nr_g*j+i] = Ge[nc_g*i+j];
230: }
231: }
232: }
234: static void BForm_Stabilisation(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
235: {
236: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
237: PetscScalar Ni_p[NODES_PER_EL],GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
238: PetscScalar J_p,fac,eta_avg;
239: PetscInt p,i,j,ngp;
241: /* define quadrature rule */
242: CreateGaussQuadrature(&ngp,gp_xi,gp_weight);
244: /* evaluate bilinear form */
245: for (p = 0; p < ngp; p++) {
246: EvaluateBasis_Q1(gp_xi[p],Ni_p);
247: EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
248: EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
249: fac = gp_weight[p]*J_p;
251: for (i = 0; i < NODES_PER_EL; i++) {
252: for (j = 0; j < NODES_PER_EL; j++) {
253: Ke[NODES_PER_EL*i+j] -= fac*(Ni_p[i]*Ni_p[j]-0.0625);
254: }
255: }
256: }
258: /* scale */
259: eta_avg = 0.0;
260: for (p = 0; p < ngp; p++) eta_avg += eta[p];
261: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
262: fac = 1.0/eta_avg;
263: for (i = 0; i < NODES_PER_EL; i++) {
264: for (j = 0; j < NODES_PER_EL; j++) {
265: Ke[NODES_PER_EL*i+j] = fac*Ke[NODES_PER_EL*i+j];
266: }
267: }
268: }
270: static void BForm_ScaledMassMatrix(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])
271: {
272: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
273: PetscScalar Ni_p[NODES_PER_EL],GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
274: PetscScalar J_p,fac,eta_avg;
275: PetscInt p,i,j,ngp;
277: /* define quadrature rule */
278: CreateGaussQuadrature(&ngp,gp_xi,gp_weight);
280: /* evaluate bilinear form */
281: for (p = 0; p < ngp; p++) {
282: EvaluateBasis_Q1(gp_xi[p],Ni_p);
283: EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
284: EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
285: fac = gp_weight[p]*J_p;
287: for (i = 0; i < NODES_PER_EL; i++) {
288: for (j = 0; j < NODES_PER_EL; j++) {
289: Ke[NODES_PER_EL*i+j] -= fac*Ni_p[i]*Ni_p[j];
290: }
291: }
292: }
294: /* scale */
295: eta_avg = 0.0;
296: for (p = 0; p < ngp; p++) eta_avg += eta[p];
297: eta_avg = (1.0/((PetscScalar)ngp))*eta_avg;
298: fac = 1.0/eta_avg;
299: for (i = 0; i < NODES_PER_EL; i++) {
300: for (j = 0; j < NODES_PER_EL; j++) {
301: Ke[NODES_PER_EL*i+j] *= fac;
302: }
303: }
304: }
306: static void LForm_MomentumRHS(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])
307: {
308: PetscScalar gp_xi[GAUSS_POINTS][NSD],gp_weight[GAUSS_POINTS];
309: PetscScalar Ni_p[NODES_PER_EL],GNi_p[NSD][NODES_PER_EL],GNx_p[NSD][NODES_PER_EL];
310: PetscScalar J_p,fac;
311: PetscInt p,i,ngp;
313: /* define quadrature rule */
314: CreateGaussQuadrature(&ngp,gp_xi,gp_weight);
316: /* evaluate linear form */
317: for (p = 0; p < ngp; p++) {
318: EvaluateBasis_Q1(gp_xi[p],Ni_p);
319: EvaluateBasisDerivatives_Q1(gp_xi[p],GNi_p);
320: EvaluateDerivatives(GNi_p,GNx_p,coords,&J_p);
321: fac = gp_weight[p]*J_p;
323: for (i = 0; i < NODES_PER_EL; i++) {
324: Fe[NSD*i] = 0.0;
325: Fe[NSD*i+1] -= fac*Ni_p[i]*fy[p];
326: }
327: }
328: }
330: static PetscErrorCode GetElementCoords(const PetscScalar _coords[],const PetscInt e2n[],PetscScalar el_coords[])
331: {
332: PetscInt i,d;
334: /* get coords for the element */
335: for (i=0; i<4; i++) {
336: for (d=0; d<NSD; d++) {
337: el_coords[NSD*i+d] = _coords[NSD * e2n[i] + d];
338: }
339: }
340: return(0);
341: }
343: static PetscErrorCode AssembleStokes_A(Mat A,DM stokes_da,DM quadrature)
344: {
345: DM cda;
346: Vec coords;
347: const PetscScalar *_coords;
348: PetscInt u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
349: PetscInt p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
350: PetscInt nel,npe,eidx;
351: const PetscInt *element_list;
352: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
353: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
354: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
355: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
356: PetscScalar el_coords[NODES_PER_EL*NSD];
357: PetscScalar *q_eta,*prop_eta;
358: PetscErrorCode ierr;
361: MatZeroEntries(A);
362: /* setup for coords */
363: DMGetCoordinateDM(stokes_da,&cda);
364: DMGetCoordinatesLocal(stokes_da,&coords);
365: VecGetArrayRead(coords,&_coords);
367: /* setup for coefficients */
368: DMSwarmGetField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
370: DMDAGetElements(stokes_da,&nel,&npe,&element_list);
371: for (eidx = 0; eidx < nel; eidx++) {
372: const PetscInt *element = &element_list[npe*eidx];
374: /* get coords for the element */
375: GetElementCoords(_coords,element,el_coords);
377: /* get coefficients for the element */
378: prop_eta = &q_eta[GAUSS_POINTS * eidx];
380: /* initialise element stiffness matrix */
381: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
382: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
383: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
384: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
386: /* form element stiffness matrix */
387: BForm_DivT(Ae,el_coords,prop_eta);
388: BForm_Grad(Ge,el_coords);
389: BForm_Div(De,el_coords);
390: BForm_Stabilisation(Ce,el_coords,prop_eta);
392: /* insert element matrix into global matrix */
393: DMDAGetElementEqnums_up(element,u_eqn,p_eqn);
394: MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
395: MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
396: MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
397: MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
398: }
399: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
400: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
402: DMSwarmRestoreField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
403: VecRestoreArrayRead(coords,&_coords);
404: return(0);
405: }
407: static PetscErrorCode AssembleStokes_PC(Mat A,DM stokes_da,DM quadrature)
408: {
409: DM cda;
410: Vec coords;
411: const PetscScalar *_coords;
412: PetscInt u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
413: PetscInt p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
414: PetscInt nel,npe,eidx;
415: const PetscInt *element_list;
416: PetscScalar Ae[NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS];
417: PetscScalar Ge[NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS];
418: PetscScalar De[NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS];
419: PetscScalar Ce[NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS];
420: PetscScalar el_coords[NODES_PER_EL*NSD];
421: PetscScalar *q_eta,*prop_eta;
422: PetscErrorCode ierr;
425: MatZeroEntries(A);
426: /* setup for coords */
427: DMGetCoordinateDM(stokes_da,&cda);
428: DMGetCoordinatesLocal(stokes_da,&coords);
429: VecGetArrayRead(coords,&_coords);
431: /* setup for coefficients */
432: DMSwarmGetField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
434: DMDAGetElements(stokes_da,&nel,&npe,&element_list);
435: for (eidx = 0; eidx < nel; eidx++) {
436: const PetscInt *element = &element_list[npe*eidx];
438: /* get coords for the element */
439: GetElementCoords(_coords,element,el_coords);
441: /* get coefficients for the element */
442: prop_eta = &q_eta[GAUSS_POINTS * eidx];
444: /* initialise element stiffness matrix */
445: PetscMemzero(Ae,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*U_DOFS);
446: PetscMemzero(Ge,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS*NODES_PER_EL*P_DOFS);
447: PetscMemzero(De,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*U_DOFS);
448: PetscMemzero(Ce,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS*NODES_PER_EL*P_DOFS);
450: /* form element stiffness matrix */
451: BForm_DivT(Ae,el_coords,prop_eta);
452: BForm_Grad(Ge,el_coords);
453: BForm_ScaledMassMatrix(Ce,el_coords,prop_eta);
455: /* insert element matrix into global matrix */
456: DMDAGetElementEqnums_up(element,u_eqn,p_eqn);
457: MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);
458: MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);
459: MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);
460: MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);
461: }
462: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
463: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
465: DMSwarmRestoreField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
466: VecRestoreArrayRead(coords,&_coords);
468: return(0);
469: }
471: static PetscErrorCode AssembleStokes_RHS(Vec F,DM stokes_da,DM quadrature)
472: {
473: DM cda;
474: Vec coords;
475: const PetscScalar *_coords;
476: PetscInt u_eqn[NODES_PER_EL*U_DOFS]; /* 2 degrees of freedom */
477: PetscInt p_eqn[NODES_PER_EL*P_DOFS]; /* 1 degrees of freedom */
478: PetscInt nel,npe,eidx,i;
479: const PetscInt *element_list;
480: PetscScalar Fe[NODES_PER_EL*U_DOFS];
481: PetscScalar He[NODES_PER_EL*P_DOFS];
482: PetscScalar el_coords[NODES_PER_EL*NSD];
483: PetscScalar *q_rhs,*prop_fy;
484: Vec local_F;
485: PetscScalar *LA_F;
486: PetscErrorCode ierr;
489: VecZeroEntries(F);
490: /* setup for coords */
491: DMGetCoordinateDM(stokes_da,&cda);
492: DMGetCoordinatesLocal(stokes_da,&coords);
493: VecGetArrayRead(coords,&_coords);
495: /* setup for coefficients */
496: DMSwarmGetField(quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
498: /* get acces to the vector */
499: DMGetLocalVector(stokes_da,&local_F);
500: VecZeroEntries(local_F);
501: VecGetArray(local_F,&LA_F);
503: DMDAGetElements(stokes_da,&nel,&npe,&element_list);
504: for (eidx = 0; eidx < nel; eidx++) {
505: const PetscInt *element = &element_list[npe*eidx];
507: /* get coords for the element */
508: GetElementCoords(_coords,element,el_coords);
510: /* get coefficients for the element */
511: prop_fy = &q_rhs[GAUSS_POINTS * eidx];
513: /* initialise element stiffness matrix */
514: PetscMemzero(Fe,sizeof(PetscScalar)*NODES_PER_EL*U_DOFS);
515: PetscMemzero(He,sizeof(PetscScalar)*NODES_PER_EL*P_DOFS);
517: /* form element stiffness matrix */
518: LForm_MomentumRHS(Fe,el_coords,NULL,prop_fy);
520: /* insert element matrix into global matrix */
521: DMDAGetElementEqnums_up(element,u_eqn,p_eqn);
523: for (i=0; i<NODES_PER_EL*U_DOFS; i++) {
524: LA_F[ u_eqn[i] ] += Fe[i];
525: }
526: }
527: DMSwarmRestoreField(quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
528: VecRestoreArrayRead(coords,&_coords);
530: VecRestoreArray(local_F,&LA_F);
531: DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);
532: DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);
533: DMRestoreLocalVector(stokes_da,&local_F);
535: return(0);
536: }
538: PetscErrorCode DMSwarmPICInsertPointsCellwise(DM dm,DM dmc,PetscInt e,PetscInt npoints,PetscReal xi[],PetscBool proximity_initialization)
539: {
540: PetscErrorCode ierr;
541: PetscInt dim,nel,npe,q,k,d,ncurr;
542: const PetscInt *element_list;
543: Vec coor;
544: const PetscScalar *_coor;
545: PetscReal **basis,*elcoor,*xp;
546: PetscReal *swarm_coor;
547: PetscInt *swarm_cellid;
550: DMGetDimension(dm,&dim);
551: DMDAGetElements(dmc,&nel,&npe,&element_list);
553: PetscMalloc1(dim*npoints,&xp);
554: PetscMalloc1(dim*npe,&elcoor);
555: PetscMalloc1(npoints,&basis);
556: for (q=0; q<npoints; q++) {
557: PetscMalloc1(npe,&basis[q]);
559: switch (dim) {
560: case 1:
561: basis[q][0] = 0.5*(1.0 - xi[dim*q+0]);
562: basis[q][1] = 0.5*(1.0 + xi[dim*q+0]);
563: break;
564: case 2:
565: basis[q][0] = 0.25*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1]);
566: basis[q][1] = 0.25*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1]);
567: basis[q][2] = 0.25*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1]);
568: basis[q][3] = 0.25*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1]);
569: break;
571: case 3:
572: basis[q][0] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
573: basis[q][1] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 - xi[dim*q+2]);
574: basis[q][2] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
575: basis[q][3] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 - xi[dim*q+2]);
576: basis[q][4] = 0.125*(1.0 - xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
577: basis[q][5] = 0.125*(1.0 + xi[dim*q+0])*(1.0 - xi[dim*q+1])*(1.0 + xi[dim*q+2]);
578: basis[q][6] = 0.125*(1.0 + xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
579: basis[q][7] = 0.125*(1.0 - xi[dim*q+0])*(1.0 + xi[dim*q+1])*(1.0 + xi[dim*q+2]);
580: break;
581: }
582: }
584: DMGetCoordinatesLocal(dmc,&coor);
585: VecGetArrayRead(coor,&_coor);
586: /* compute and store the coordinates for the new points */
587: {
588: const PetscInt *element = &element_list[npe*e];
590: for (k=0; k<npe; k++) {
591: for (d=0; d<dim; d++) {
592: elcoor[dim*k+d] = PetscRealPart(_coor[ dim*element[k] + d ]);
593: }
594: }
595: for (q=0; q<npoints; q++) {
596: for (d=0; d<dim; d++) {
597: xp[dim*q+d] = 0.0;
598: }
599: for (k=0; k<npe; k++) {
600: for (d=0; d<dim; d++) {
601: xp[dim*q+d] += basis[q][k] * elcoor[dim*k+d];
602: }
603: }
604: }
605: }
606: VecRestoreArrayRead(coor,&_coor);
607: DMDARestoreElements(dmc,&nel,&npe,&element_list);
609: DMSwarmGetLocalSize(dm,&ncurr);
610: DMSwarmAddNPoints(dm,npoints);
612: if (proximity_initialization) {
613: PetscInt *nnlist;
614: PetscReal *coor_q,*coor_qn;
615: PetscInt npoints_e,*plist_e;
617: DMSwarmSortGetPointsPerCell(dm,e,&npoints_e,&plist_e);
619: PetscMalloc1(npoints,&nnlist);
620: /* find nearest neighour points in this cell */
621: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
622: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
623: for (q=0; q<npoints; q++) {
624: PetscInt qn,nearest_neighour = -1;
625: PetscReal sep,min_sep = PETSC_MAX_REAL;
627: coor_q = &xp[dim*q];
628: for (qn=0; qn<npoints_e; qn++) {
629: coor_qn = &swarm_coor[dim*plist_e[qn]];
630: sep = 0.0;
631: for (d=0; d<dim; d++) {
632: sep += (coor_q[d]-coor_qn[d])*(coor_q[d]-coor_qn[d]);
633: }
634: if (sep < min_sep) {
635: nearest_neighour = plist_e[qn];
636: min_sep = sep;
637: }
638: }
639: if (nearest_neighour == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cell %D is empty - cannot initalize using nearest neighbours",e);
640: nnlist[q] = nearest_neighour;
641: }
642: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
643: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
645: /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
646: for (q=0; q<npoints; q++) {
647: DMSwarmCopyPoint(dm,nnlist[q],ncurr+q);
648: }
649: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
650: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
651: for (q=0; q<npoints; q++) {
652: /* set the coordinates */
653: for (d=0; d<dim; d++) {
654: swarm_coor[dim*(ncurr+q)+d] = xp[dim*q+d];
655: }
656: /* set the cell index */
657: swarm_cellid[ncurr+q] = e;
658: }
659: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
660: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
662: PetscFree(plist_e);
663: PetscFree(nnlist);
664: } else {
665: DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
666: DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
667: for (q=0; q<npoints; q++) {
668: /* set the coordinates */
669: for (d=0; d<dim; d++) {
670: swarm_coor[dim*(ncurr+q)+d] = xp[dim*q+d];
671: }
672: /* set the cell index */
673: swarm_cellid[ncurr+q] = e;
674: }
675: DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);
676: DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);
677: }
679: PetscFree(xp);
680: PetscFree(elcoor);
681: for (q=0; q<npoints; q++) {
682: PetscFree(basis[q]);
683: }
684: PetscFree(basis);
685: return(0);
686: }
688: PetscErrorCode MaterialPoint_PopulateCell(DM dm_vp,DM dm_mpoint)
689: {
690: PetscInt _npe,_nel,e,nel;
691: const PetscInt *element;
692: DM dmc;
693: PetscQuadrature quadrature;
694: const PetscReal *xi;
695: PetscInt npoints_q,cnt,cnt_g;
696: PetscErrorCode ierr;
699: DMDAGetElements(dm_vp,&_nel,&_npe,&element);
700: nel = _nel;
701: DMDARestoreElements(dm_vp,&_nel,&_npe,&element);
703: PetscDTGaussTensorQuadrature(2,1,4,-1.0,1.0,&quadrature);
704: PetscQuadratureGetData(quadrature,NULL,NULL,&npoints_q,&xi,NULL);
705: DMSwarmGetCellDM(dm_mpoint,&dmc);
707: DMSwarmSortGetAccess(dm_mpoint);
709: cnt = 0;
710: for (e=0; e<nel; e++) {
711: PetscInt npoints_per_cell;
713: DMSwarmSortGetNumberOfPointsPerCell(dm_mpoint,e,&npoints_per_cell);
715: if (npoints_per_cell < 12) {
716: DMSwarmPICInsertPointsCellwise(dm_mpoint,dm_vp,e,npoints_q,(PetscReal*)xi,PETSC_TRUE);
717: cnt++;
718: }
719: }
720: MPI_Allreduce(&cnt,&cnt_g,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
721: if (cnt_g > 0) {
722: PetscPrintf(PETSC_COMM_WORLD,".... ....pop cont: adjusted %D cells\n",cnt_g);
723: }
725: DMSwarmSortRestoreAccess(dm_mpoint);
726: PetscQuadratureDestroy(&quadrature);
727: return(0);
728: }
730: PetscErrorCode MaterialPoint_AdvectRK1(DM dm_vp,Vec vp,PetscReal dt,DM dm_mpoint)
731: {
732: PetscErrorCode ierr;
733: Vec vp_l,coor_l;
734: const PetscScalar *LA_vp;
735: PetscInt i,p,e,npoints,nel,npe;
736: PetscInt *mpfield_cell;
737: PetscReal *mpfield_coor;
738: const PetscInt *element_list;
739: const PetscInt *element;
740: PetscScalar xi_p[NSD],Ni[NODES_PER_EL];
741: const PetscScalar *LA_coor;
742: PetscScalar dx[NSD];
745: DMGetCoordinatesLocal(dm_vp,&coor_l);
746: VecGetArrayRead(coor_l,&LA_coor);
748: DMGetLocalVector(dm_vp,&vp_l);
749: DMGlobalToLocalBegin(dm_vp,vp,INSERT_VALUES,vp_l);
750: DMGlobalToLocalEnd(dm_vp,vp,INSERT_VALUES,vp_l);
751: VecGetArrayRead(vp_l,&LA_vp);
753: DMDAGetElements(dm_vp,&nel,&npe,&element_list);
754: DMSwarmGetLocalSize(dm_mpoint,&npoints);
755: DMSwarmGetField(dm_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
756: DMSwarmGetField(dm_mpoint,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
757: for (p=0; p<npoints; p++) {
758: PetscReal *coor_p;
759: PetscScalar vel_n[NSD*NODES_PER_EL],vel_p[NSD];
760: const PetscScalar *x0;
761: const PetscScalar *x2;
763: e = mpfield_cell[p];
764: coor_p = &mpfield_coor[NSD*p];
765: element = &element_list[NODES_PER_EL*e];
767: /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
768: x0 = &LA_coor[NSD*element[0]];
769: x2 = &LA_coor[NSD*element[2]];
771: dx[0] = x2[0] - x0[0];
772: dx[1] = x2[1] - x0[1];
774: xi_p[0] = 2.0 * (coor_p[0] - x0[0])/dx[0] - 1.0;
775: xi_p[1] = 2.0 * (coor_p[1] - x0[1])/dx[1] - 1.0;
776: 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);
777: 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);
778: 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);
779: 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);
781: /* evaluate basis functions */
782: EvaluateBasis_Q1(xi_p,Ni);
784: /* get cell nodal velocities */
785: for (i=0; i<NODES_PER_EL; i++) {
786: PetscInt nid;
788: nid = element[i];
789: vel_n[NSD*i+0] = LA_vp[(NSD+1)*nid+0];
790: vel_n[NSD*i+1] = LA_vp[(NSD+1)*nid+1];
791: }
793: /* interpolate velocity */
794: vel_p[0] = vel_p[1] = 0.0;
795: for (i=0; i<NODES_PER_EL; i++) {
796: vel_p[0] += Ni[i] * vel_n[NSD*i+0];
797: vel_p[1] += Ni[i] * vel_n[NSD*i+1];
798: }
800: coor_p[0] += dt * PetscRealPart(vel_p[0]);
801: coor_p[1] += dt * PetscRealPart(vel_p[1]);
802: }
804: DMSwarmRestoreField(dm_mpoint,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);
805: DMSwarmRestoreField(dm_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);
806: DMDARestoreElements(dm_vp,&nel,&npe,&element_list);
807: VecRestoreArrayRead(vp_l,&LA_vp);
808: DMRestoreLocalVector(dm_vp,&vp_l);
809: VecRestoreArrayRead(coor_l,&LA_coor);
810: return(0);
811: }
813: PetscErrorCode MaterialPoint_Interpolate(DM dm,Vec eta_v,Vec rho_v,DM dm_quadrature)
814: {
815: Vec eta_l,rho_l;
816: PetscScalar *_eta_l,*_rho_l;
817: PetscInt nqp,npe,nel;
818: PetscScalar qp_xi[GAUSS_POINTS][NSD];
819: PetscScalar qp_weight[GAUSS_POINTS];
820: PetscInt q,k,e;
821: PetscScalar Ni[GAUSS_POINTS][NODES_PER_EL];
822: const PetscInt *element_list;
823: PetscReal *q_eta,*q_rhs;
827: /* define quadrature rule */
828: CreateGaussQuadrature(&nqp,qp_xi,qp_weight);
829: for (q=0; q<nqp; q++) {
830: EvaluateBasis_Q1(qp_xi[q],Ni[q]);
831: }
833: DMGetLocalVector(dm,&eta_l);
834: DMGetLocalVector(dm,&rho_l);
836: DMGlobalToLocalBegin(dm,eta_v,INSERT_VALUES,eta_l);
837: DMGlobalToLocalEnd(dm,eta_v,INSERT_VALUES,eta_l);
838: DMGlobalToLocalBegin(dm,rho_v,INSERT_VALUES,rho_l);
839: DMGlobalToLocalEnd(dm,rho_v,INSERT_VALUES,rho_l);
841: VecGetArray(eta_l,&_eta_l);
842: VecGetArray(rho_l,&_rho_l);
844: DMSwarmGetField(dm_quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
845: DMSwarmGetField(dm_quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
847: DMDAGetElements(dm,&nel,&npe,&element_list);
848: for (e=0; e<nel; e++) {
849: PetscScalar eta_field_e[NODES_PER_EL];
850: PetscScalar rho_field_e[NODES_PER_EL];
851: const PetscInt *element = &element_list[4*e];
853: for (k=0; k<NODES_PER_EL; k++) {
854: eta_field_e[k] = _eta_l[ element[k] ];
855: rho_field_e[k] = _rho_l[ element[k] ];
856: }
858: for (q=0; q<nqp; q++) {
859: PetscScalar eta_q,rho_q;
861: eta_q = rho_q = 0.0;
862: for (k=0; k<NODES_PER_EL; k++) {
863: eta_q += Ni[q][k] * eta_field_e[k];
864: rho_q += Ni[q][k] * rho_field_e[k];
865: }
867: q_eta[nqp*e+q] = PetscRealPart(eta_q);
868: q_rhs[nqp*e+q] = PetscRealPart(rho_q);
869: }
870: }
871: DMDARestoreElements(dm,&nel,&npe,&element_list);
873: DMSwarmRestoreField(dm_quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);
874: DMSwarmRestoreField(dm_quadrature,"eta_q",NULL,NULL,(void**)&q_eta);
876: VecRestoreArray(rho_l,&_rho_l);
877: VecRestoreArray(eta_l,&_eta_l);
878: DMRestoreLocalVector(dm,&rho_l);
879: DMRestoreLocalVector(dm,&eta_l);
880: return(0);
881: }
883: static PetscErrorCode SolveTimeDepStokes(PetscInt mx,PetscInt my)
884: {
885: DM dm_stokes,dm_coeff;
886: PetscInt u_dof,p_dof,dof,stencil_width;
887: Mat A,B;
888: PetscInt nel_local;
889: Vec eta_v,rho_v;
890: Vec f,X;
891: KSP ksp;
892: PC pc;
893: char filename[PETSC_MAX_PATH_LEN];
894: DM dms_quadrature,dms_mpoint;
895: PetscInt nel,npe,npoints;
896: const PetscInt *element_list;
897: PetscInt tk,nt,dump_freq;
898: PetscReal dt,dt_max = 0.0;
899: PetscReal vx[2],vy[2],max_v = 0.0,max_v_step,dh;
900: PetscErrorCode ierr;
901: const char *fieldnames[] = { "eta" , "rho" };
902: Vec *pfields;
903: PetscInt ppcell = 1;
904: PetscReal time,delta_eta = 1.0;
905: PetscBool randomize_coords = PETSC_FALSE;
906: PetscReal randomize_fac = 0.25;
907: PetscBool no_view = PETSC_FALSE;
908: PetscBool isbddc;
911: /*
912: Generate the DMDA for the velocity and pressure spaces.
913: We use Q1 elements for both fields.
914: The Q1 FE basis on a regular mesh has a 9-point stencil (DMDA_STENCIL_BOX)
915: The number of nodes in each direction is mx+1, my+1
916: */
917: u_dof = U_DOFS; /* Vx, Vy - velocities */
918: p_dof = P_DOFS; /* p - pressure */
919: dof = u_dof + p_dof;
920: stencil_width = 1;
921: 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);
922: DMDASetElementType(dm_stokes,DMDA_ELEMENT_Q1);
923: DMSetMatType(dm_stokes,MATAIJ);
924: DMSetFromOptions(dm_stokes);
925: DMSetUp(dm_stokes);
926: DMDASetFieldName(dm_stokes,0,"ux");
927: DMDASetFieldName(dm_stokes,1,"uy");
928: DMDASetFieldName(dm_stokes,2,"p");
930: /* unit box [0,0.9142] x [0,1] */
931: DMDASetUniformCoordinates(dm_stokes,0.0,0.9142,0.0,1.0,0.,0.);
932: dh = 1.0/((PetscReal)(mx));
934: /* Get local number of elements */
935: {
936: DMDAGetElements(dm_stokes,&nel,&npe,&element_list);
938: nel_local = nel;
940: DMDARestoreElements(dm_stokes,&nel,&npe,&element_list);
941: }
943: /* Create DMDA for representing scalar fields */
944: DMDACreateCompatibleDMDA(dm_stokes,1,&dm_coeff);
946: /* Create the swarm for storing quadrature point values */
947: DMCreate(PETSC_COMM_WORLD,&dms_quadrature);
948: DMSetType(dms_quadrature,DMSWARM);
949: DMSetDimension(dms_quadrature,2);
951: /* Register fields for viscosity and density on the quadrature points */
952: DMSwarmRegisterPetscDatatypeField(dms_quadrature,"eta_q",1,PETSC_REAL);
953: DMSwarmRegisterPetscDatatypeField(dms_quadrature,"rho_q",1,PETSC_REAL);
954: DMSwarmFinalizeFieldRegister(dms_quadrature);
955: DMSwarmSetLocalSizes(dms_quadrature,nel_local * GAUSS_POINTS,0);
957: /* Create the material point swarm */
958: DMCreate(PETSC_COMM_WORLD,&dms_mpoint);
959: DMSetType(dms_mpoint,DMSWARM);
960: DMSetDimension(dms_mpoint,2);
962: /* Configure the material point swarm to be of type Particle-In-Cell */
963: DMSwarmSetType(dms_mpoint,DMSWARM_PIC);
965: /*
966: Specify the DM to use for point location and projections
967: within the context of a PIC scheme
968: */
969: DMSwarmSetCellDM(dms_mpoint,dm_coeff);
971: /* Register fields for viscosity and density */
972: DMSwarmRegisterPetscDatatypeField(dms_mpoint,"eta",1,PETSC_REAL);
973: DMSwarmRegisterPetscDatatypeField(dms_mpoint,"rho",1,PETSC_REAL);
974: DMSwarmFinalizeFieldRegister(dms_mpoint);
976: PetscOptionsGetInt(NULL,NULL,"-ppcell",&ppcell,NULL);
977: DMSwarmSetLocalSizes(dms_mpoint,nel_local * ppcell,100);
979: /*
980: Layout the material points in space using the cell DM.
981: Particle coordinates are defined by cell wise using different methods.
982: - DMSWARMPIC_LAYOUT_GAUSS defines particles coordinates at the positions
983: corresponding to a Gauss quadrature rule with
984: ppcell points in each direction.
985: - DMSWARMPIC_LAYOUT_REGULAR defines particle coordinates at the centoid of
986: ppcell x ppcell quadralaterals defined within the
987: reference element.
988: - DMSWARMPIC_LAYOUT_SUBDIVISION defines particles coordinates at the centroid
989: of each quadralateral obtained by sub-dividing
990: the reference element cell ppcell times.
991: */
992: DMSwarmInsertPointsUsingCellDM(dms_mpoint,DMSWARMPIC_LAYOUT_SUBDIVISION,ppcell);
994: /*
995: Defne a high resolution layer of material points across the material interface
996: */
997: {
998: PetscInt npoints_dir_x[2];
999: PetscReal min[2],max[2];
1001: npoints_dir_x[0] = (PetscInt)(0.9142/(0.05*dh));
1002: npoints_dir_x[1] = (PetscInt)((0.25-0.15)/(0.05*dh));
1003: min[0] = 0.0; max[0] = 0.9142;
1004: min[1] = 0.05; max[1] = 0.35;
1005: DMSwarmSetPointsUniformCoordinates(dms_mpoint,min,max,npoints_dir_x,ADD_VALUES);
1006: }
1008: /*
1009: Define a high resolution layer of material points near the surface of the domain
1010: to deal with weakly compressible Q1-Q1 elements. These elements "self compact"
1011: when applied to buouyancy driven flow. The error in div(u) is O(h).
1012: */
1013: {
1014: PetscInt npoints_dir_x[2];
1015: PetscReal min[2],max[2];
1017: npoints_dir_x[0] = (PetscInt)(0.9142/(0.25*dh));
1018: npoints_dir_x[1] = (PetscInt)(3.0*dh/(0.25*dh));
1019: min[0] = 0.0; max[0] = 0.9142;
1020: min[1] = 1.0 - 3.0*dh; max[1] = 1.0-0.0001;
1021: DMSwarmSetPointsUniformCoordinates(dms_mpoint,min,max,npoints_dir_x,ADD_VALUES);
1022: }
1024: DMView(dms_mpoint,PETSC_VIEWER_STDOUT_WORLD);
1026: /* Define initial material properties on each particle in the material point swarm */
1027: PetscOptionsGetReal(NULL,NULL,"-delta_eta",&delta_eta,NULL);
1028: PetscOptionsGetBool(NULL,NULL,"-randomize_coords",&randomize_coords,NULL);
1029: PetscOptionsGetReal(NULL,NULL,"-randomize_fac",&randomize_fac,NULL);
1030: if (randomize_fac > 1.0) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"The value of -randomize_fac should be <= 1.0");
1031: {
1032: PetscReal *array_x,*array_e,*array_r;
1033: PetscInt p;
1034: PetscRandom r;
1035: PetscMPIInt rank;
1037: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1039: PetscRandomCreate(PETSC_COMM_SELF,&r);
1040: PetscRandomSetInterval(r,-randomize_fac*dh,randomize_fac*dh);
1041: PetscRandomSetSeed(r,(unsigned long)rank);
1042: PetscRandomSeed(r);
1044: DMDAGetElements(dm_stokes,&nel,&npe,&element_list);
1046: /*
1047: Fetch the registered data from the material point DMSwarm.
1048: The fields "eta" and "rho" were registered by this example.
1049: The field identified by the the variable DMSwarmPICField_coor
1050: was registered by the DMSwarm implementation when the function
1051: DMSwarmSetType(dms_mpoint,DMSWARM_PIC)
1052: was called. The returned array defines the coordinates of each
1053: material point in the point swarm.
1054: */
1055: DMSwarmGetField(dms_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&array_x);
1056: DMSwarmGetField(dms_mpoint,"eta", NULL,NULL,(void**)&array_e);
1057: DMSwarmGetField(dms_mpoint,"rho", NULL,NULL,(void**)&array_r);
1059: DMSwarmGetLocalSize(dms_mpoint,&npoints);
1060: for (p = 0; p < npoints; p++) {
1061: PetscReal x_p[2],rr[2];
1063: if (randomize_coords) {
1064: PetscRandomGetValueReal(r,&rr[0]);
1065: PetscRandomGetValueReal(r,&rr[1]);
1066: array_x[2*p + 0] += rr[0];
1067: array_x[2*p + 1] += rr[1];
1068: }
1070: /* Get the coordinates of point, p */
1071: x_p[0] = array_x[2*p + 0];
1072: x_p[1] = array_x[2*p + 1];
1074: if (x_p[1] < (0.2 + 0.02*PetscCosReal(PETSC_PI*x_p[0]/0.9142))) {
1075: /* Material properties below the interface */
1076: array_e[p] = 1.0 * (1.0/delta_eta);
1077: array_r[p] = 0.0;
1078: } else {
1079: /* Material properties above the interface */
1080: array_e[p] = 1.0;
1081: array_r[p] = 1.0;
1082: }
1083: }
1085: /*
1086: Restore the fetched data fields from the material point DMSwarm.
1087: Calling the Restore function invalidates the points array_r, array_e, array_x
1088: by setting them to NULL.
1089: */
1090: DMSwarmRestoreField(dms_mpoint,"rho",NULL,NULL,(void**)&array_r);
1091: DMSwarmRestoreField(dms_mpoint,"eta",NULL,NULL,(void**)&array_e);
1092: DMSwarmRestoreField(dms_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&array_x);
1094: DMDARestoreElements(dm_stokes,&nel,&npe,&element_list);
1095: PetscRandomDestroy(&r);
1096: }
1098: /*
1099: If the particle coordinates where randomly shifted, they may have crossed into another
1100: element, or into another sub-domain. To account for this we call the Migrate function.
1101: */
1102: if (randomize_coords) {
1103: DMSwarmMigrate(dms_mpoint,PETSC_TRUE);
1104: }
1106: PetscOptionsGetBool(NULL,NULL,"-no_view",&no_view,NULL);
1107: if (!no_view) {
1108: DMSwarmViewXDMF(dms_mpoint,"ic_coeff_dms.xmf");
1109: }
1111: /* project the swarm properties */
1112: DMSwarmProjectFields(dms_mpoint,2,fieldnames,&pfields,PETSC_FALSE);
1113: eta_v = pfields[0];
1114: rho_v = pfields[1];
1115: PetscObjectSetName((PetscObject)eta_v,"eta");
1116: PetscObjectSetName((PetscObject)rho_v,"rho");
1117: MaterialPoint_Interpolate(dm_coeff,eta_v,rho_v,dms_quadrature);
1119: /* view projected coefficients eta and rho */
1120: if (!no_view) {
1121: PetscViewer viewer;
1123: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1124: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1125: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1126: PetscViewerFileSetName(viewer,"ic_coeff_dmda.vts");
1127: VecView(eta_v,viewer);
1128: VecView(rho_v,viewer);
1129: PetscViewerDestroy(&viewer);
1130: }
1132: DMCreateMatrix(dm_stokes,&A);
1133: DMCreateMatrix(dm_stokes,&B);
1134: DMCreateGlobalVector(dm_stokes,&f);
1135: DMCreateGlobalVector(dm_stokes,&X);
1137: AssembleStokes_A(A,dm_stokes,dms_quadrature);
1138: AssembleStokes_PC(B,dm_stokes,dms_quadrature);
1139: AssembleStokes_RHS(f,dm_stokes,dms_quadrature);
1141: DMDAApplyBoundaryConditions(dm_stokes,A,f);
1142: DMDAApplyBoundaryConditions(dm_stokes,B,NULL);
1144: KSPCreate(PETSC_COMM_WORLD,&ksp);
1145: KSPSetOptionsPrefix(ksp,"stokes_");
1146: KSPSetDM(ksp,dm_stokes);
1147: KSPSetDMActive(ksp,PETSC_FALSE);
1148: KSPSetOperators(ksp,A,B);
1149: KSPSetFromOptions(ksp);
1150: KSPGetPC(ksp,&pc);
1151: PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);
1152: if (isbddc) {
1153: KSPSetOperators(ksp,A,A);
1154: }
1156: /* Define u-v-p indices for fieldsplit */
1157: {
1158: PC pc;
1159: const PetscInt ufields[] = {0,1},pfields[1] = {2};
1161: KSPGetPC(ksp,&pc);
1162: PCFieldSplitSetBlockSize(pc,3);
1163: PCFieldSplitSetFields(pc,"u",2,ufields,ufields);
1164: PCFieldSplitSetFields(pc,"p",1,pfields,pfields);
1165: }
1167: /* If using a fieldsplit preconditioner, attach a DMDA to the velocity split so that geometric multigrid can be used */
1168: {
1169: PC pc,pc_u;
1170: KSP *sub_ksp,ksp_u;
1171: PetscInt nsplits;
1172: DM dm_u;
1173: PetscBool is_pcfs;
1175: KSPGetPC(ksp,&pc);
1177: is_pcfs = PETSC_FALSE;
1178: PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&is_pcfs);
1180: if (is_pcfs) {
1181: KSPSetUp(ksp);
1182: KSPGetPC(ksp,&pc);
1183: PCFieldSplitGetSubKSP(pc,&nsplits,&sub_ksp);
1184: ksp_u = sub_ksp[0];
1185: PetscFree(sub_ksp);
1187: if (nsplits == 2) {
1188: DMDACreateCompatibleDMDA(dm_stokes,2,&dm_u);
1190: KSPSetDM(ksp_u,dm_u);
1191: KSPSetDMActive(ksp_u,PETSC_FALSE);
1192: DMDestroy(&dm_u);
1194: /* enforce galerkin coarse grids be used */
1195: KSPGetPC(ksp_u,&pc_u);
1196: PCMGSetGalerkin(pc_u,PC_MG_GALERKIN_PMAT);
1197: }
1198: }
1199: }
1201: dump_freq = 10;
1202: PetscOptionsGetInt(NULL,NULL,"-dump_freq",&dump_freq,NULL);
1203: nt = 10;
1204: PetscOptionsGetInt(NULL,NULL,"-nt",&nt,NULL);
1205: time = 0.0;
1206: for (tk=1; tk<=nt; tk++) {
1208: PetscPrintf(PETSC_COMM_WORLD,".... assemble\n");
1209: AssembleStokes_A(A,dm_stokes,dms_quadrature);
1210: AssembleStokes_PC(B,dm_stokes,dms_quadrature);
1211: AssembleStokes_RHS(f,dm_stokes,dms_quadrature);
1213: PetscPrintf(PETSC_COMM_WORLD,".... bc imposition\n");
1214: DMDAApplyBoundaryConditions(dm_stokes,A,f);
1215: DMDAApplyBoundaryConditions(dm_stokes,B,NULL);
1217: PetscPrintf(PETSC_COMM_WORLD,".... solve\n");
1218: KSPSetOperators(ksp,A, isbddc ? A : B);
1219: KSPSolve(ksp,f,X);
1221: VecStrideMax(X,0,NULL,&vx[1]);
1222: VecStrideMax(X,1,NULL,&vy[1]);
1223: VecStrideMin(X,0,NULL,&vx[0]);
1224: VecStrideMin(X,1,NULL,&vy[0]);
1226: max_v_step = PetscMax(vx[0],vx[1]);
1227: max_v_step = PetscMax(max_v_step,vy[0]);
1228: max_v_step = PetscMax(max_v_step,vy[1]);
1229: max_v = PetscMax(max_v,max_v_step);
1231: dt_max = 2.0;
1232: dt = 0.5 * (dh / max_v_step);
1233: 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);
1234: dt = PetscMin(dt_max,dt);
1236: /* advect */
1237: PetscPrintf(PETSC_COMM_WORLD,".... advect\n");
1238: MaterialPoint_AdvectRK1(dm_stokes,X,dt,dms_mpoint);
1240: /* migrate */
1241: PetscPrintf(PETSC_COMM_WORLD,".... migrate\n");
1242: DMSwarmMigrate(dms_mpoint,PETSC_TRUE);
1244: /* update cell population */
1245: PetscPrintf(PETSC_COMM_WORLD,".... populate cells\n");
1246: MaterialPoint_PopulateCell(dm_stokes,dms_mpoint);
1248: /* update coefficients on quadrature points */
1249: PetscPrintf(PETSC_COMM_WORLD,".... project\n");
1250: DMSwarmProjectFields(dms_mpoint,2,fieldnames,&pfields,PETSC_TRUE);
1251: eta_v = pfields[0];
1252: rho_v = pfields[1];
1253: PetscPrintf(PETSC_COMM_WORLD,".... interp\n");
1254: MaterialPoint_Interpolate(dm_coeff,eta_v,rho_v,dms_quadrature);
1256: if (tk%dump_freq == 0) {
1257: PetscViewer viewer;
1259: PetscPrintf(PETSC_COMM_WORLD,".... write XDMF, VTS\n");
1260: PetscSNPrintf(filename,PETSC_MAX_PATH_LEN-1,"step%.4D_coeff_dms.xmf",tk);
1261: DMSwarmViewXDMF(dms_mpoint,filename);
1263: PetscSNPrintf(filename,PETSC_MAX_PATH_LEN-1,"step%.4D_vp_dm.vts",tk);
1264: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
1265: PetscViewerSetType(viewer,PETSCVIEWERVTK);
1266: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1267: PetscViewerFileSetName(viewer,filename);
1268: VecView(X,viewer);
1269: PetscViewerDestroy(&viewer);
1270: }
1271: time += dt;
1272: PetscPrintf(PETSC_COMM_WORLD,"step %D : time %1.2e \n",tk,(double)time);
1273: }
1275: KSPDestroy(&ksp);
1276: VecDestroy(&X);
1277: VecDestroy(&f);
1278: MatDestroy(&A);
1279: MatDestroy(&B);
1280: VecDestroy(&eta_v);
1281: VecDestroy(&rho_v);
1282: PetscFree(pfields);
1284: DMDestroy(&dms_mpoint);
1285: DMDestroy(&dms_quadrature);
1286: DMDestroy(&dm_coeff);
1287: DMDestroy(&dm_stokes);
1288: return(0);
1289: }
1291: /*
1292: <sequential run>
1293: ./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
1294: */
1295: int main(int argc,char **args)
1296: {
1298: PetscInt mx,my;
1299: PetscBool set = PETSC_FALSE;
1301: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
1302: mx = my = 10;
1303: PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);
1304: PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);
1305: PetscOptionsGetInt(NULL,NULL,"-mxy",&mx,&set);
1306: if (set) {
1307: my = mx;
1308: }
1309: SolveTimeDepStokes(mx,my);
1310: PetscFinalize();
1311: return ierr;
1312: }
1314: /* -------------------------- helpers for boundary conditions -------------------------------- */
1315: static PetscErrorCode BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)
1316: {
1317: DM cda;
1318: Vec coords;
1319: PetscInt si,sj,nx,ny,i,j;
1320: PetscInt M,N;
1321: DMDACoor2d **_coords;
1322: const PetscInt *g_idx;
1323: PetscInt *bc_global_ids;
1324: PetscScalar *bc_vals;
1325: PetscInt nbcs;
1326: PetscInt n_dofs;
1327: PetscErrorCode ierr;
1328: ISLocalToGlobalMapping ltogm;
1331: DMGetLocalToGlobalMapping(da,<ogm);
1332: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1334: DMGetCoordinateDM(da,&cda);
1335: DMGetCoordinatesLocal(da,&coords);
1336: DMDAVecGetArray(cda,coords,&_coords);
1337: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1338: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1340: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1341: PetscMalloc1(ny*n_dofs,&bc_vals);
1343: /* init the entries to -1 so VecSetValues will ignore them */
1344: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1346: i = nx-1;
1347: for (j = 0; j < ny; j++) {
1348: PetscInt local_id;
1350: local_id = i+j*nx;
1352: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1354: bc_vals[j] = 0.0;
1355: }
1356: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1357: nbcs = 0;
1358: if ((si+nx) == (M)) nbcs = ny;
1360: if (b) {
1361: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1362: VecAssemblyBegin(b);
1363: VecAssemblyEnd(b);
1364: }
1365: if (A) {
1366: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1367: }
1369: PetscFree(bc_vals);
1370: PetscFree(bc_global_ids);
1372: DMDAVecRestoreArray(cda,coords,&_coords);
1373: return(0);
1374: }
1376: static PetscErrorCode BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)
1377: {
1378: DM cda;
1379: Vec coords;
1380: PetscInt si,sj,nx,ny,i,j;
1381: PetscInt M,N;
1382: DMDACoor2d **_coords;
1383: const PetscInt *g_idx;
1384: PetscInt *bc_global_ids;
1385: PetscScalar *bc_vals;
1386: PetscInt nbcs;
1387: PetscInt n_dofs;
1388: PetscErrorCode ierr;
1389: ISLocalToGlobalMapping ltogm;
1392: DMGetLocalToGlobalMapping(da,<ogm);
1393: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1395: DMGetCoordinateDM(da,&cda);
1396: DMGetCoordinatesLocal(da,&coords);
1397: DMDAVecGetArray(cda,coords,&_coords);
1398: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1399: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1401: PetscMalloc1(ny*n_dofs,&bc_global_ids);
1402: PetscMalloc1(ny*n_dofs,&bc_vals);
1404: /* init the entries to -1 so VecSetValues will ignore them */
1405: for (i = 0; i < ny*n_dofs; i++) bc_global_ids[i] = -1;
1407: i = 0;
1408: for (j = 0; j < ny; j++) {
1409: PetscInt local_id;
1411: local_id = i+j*nx;
1413: bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1415: bc_vals[j] = 0.0;
1416: }
1417: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1418: nbcs = 0;
1419: if (si == 0) nbcs = ny;
1421: if (b) {
1422: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1423: VecAssemblyBegin(b);
1424: VecAssemblyEnd(b);
1425: }
1427: if (A) {
1428: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1429: }
1431: PetscFree(bc_vals);
1432: PetscFree(bc_global_ids);
1434: DMDAVecRestoreArray(cda,coords,&_coords);
1435: return(0);
1436: }
1438: static PetscErrorCode BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)
1439: {
1440: DM cda;
1441: Vec coords;
1442: PetscInt si,sj,nx,ny,i,j;
1443: PetscInt M,N;
1444: DMDACoor2d **_coords;
1445: const PetscInt *g_idx;
1446: PetscInt *bc_global_ids;
1447: PetscScalar *bc_vals;
1448: PetscInt nbcs;
1449: PetscInt n_dofs;
1450: PetscErrorCode ierr;
1451: ISLocalToGlobalMapping ltogm;
1454: DMGetLocalToGlobalMapping(da,<ogm);
1455: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1457: DMGetCoordinateDM(da,&cda);
1458: DMGetCoordinatesLocal(da,&coords);
1459: DMDAVecGetArray(cda,coords,&_coords);
1460: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1461: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1463: PetscMalloc1(nx,&bc_global_ids);
1464: PetscMalloc1(nx,&bc_vals);
1466: /* init the entries to -1 so VecSetValues will ignore them */
1467: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1469: j = ny-1;
1470: for (i = 0; i < nx; i++) {
1471: PetscInt local_id;
1473: local_id = i+j*nx;
1475: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1477: bc_vals[i] = 0.0;
1478: }
1479: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1480: nbcs = 0;
1481: if ((sj+ny) == (N)) nbcs = nx;
1483: if (b) {
1484: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1485: VecAssemblyBegin(b);
1486: VecAssemblyEnd(b);
1487: }
1488: if (A) {
1489: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1490: }
1492: PetscFree(bc_vals);
1493: PetscFree(bc_global_ids);
1495: DMDAVecRestoreArray(cda,coords,&_coords);
1496: return(0);
1497: }
1499: static PetscErrorCode BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)
1500: {
1501: DM cda;
1502: Vec coords;
1503: PetscInt si,sj,nx,ny,i,j;
1504: PetscInt M,N;
1505: DMDACoor2d **_coords;
1506: const PetscInt *g_idx;
1507: PetscInt *bc_global_ids;
1508: PetscScalar *bc_vals;
1509: PetscInt nbcs;
1510: PetscInt n_dofs;
1511: PetscErrorCode ierr;
1512: ISLocalToGlobalMapping ltogm;
1515: DMGetLocalToGlobalMapping(da,<ogm);
1516: ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);
1518: DMGetCoordinateDM(da,&cda);
1519: DMGetCoordinatesLocal(da,&coords);
1520: DMDAVecGetArray(cda,coords,&_coords);
1521: DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);
1522: DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);
1524: PetscMalloc1(nx,&bc_global_ids);
1525: PetscMalloc1(nx,&bc_vals);
1527: /* init the entries to -1 so VecSetValues will ignore them */
1528: for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1530: j = 0;
1531: for (i = 0; i < nx; i++) {
1532: PetscInt local_id;
1534: local_id = i+j*nx;
1536: bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1538: bc_vals[i] = 0.0;
1539: }
1540: ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);
1541: nbcs = 0;
1542: if (sj == 0) nbcs = nx;
1544: if (b) {
1545: VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);
1546: VecAssemblyBegin(b);
1547: VecAssemblyEnd(b);
1548: }
1549: if (A) {
1550: MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);
1551: }
1553: PetscFree(bc_vals);
1554: PetscFree(bc_global_ids);
1556: DMDAVecRestoreArray(cda,coords,&_coords);
1557: return(0);
1558: }
1560: /*
1561: Impose free slip boundary conditions on the left/right faces: u_i n_i = 0, tau_{ij} t_j = 0
1562: Impose no slip boundray conditions on the top/bottom faces: u_i n_i = 0, u_i t_i = 0
1563: */
1564: static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes,Mat A,Vec f)
1565: {
1569: BCApplyZero_NORTH(dm_stokes,0,A,f);
1570: BCApplyZero_NORTH(dm_stokes,1,A,f);
1571: BCApplyZero_EAST(dm_stokes,0,A,f);
1572: BCApplyZero_SOUTH(dm_stokes,0,A,f);
1573: BCApplyZero_SOUTH(dm_stokes,1,A,f);
1574: BCApplyZero_WEST(dm_stokes,0,A,f);
1575: return(0);
1576: }
1578: /*TEST
1580: test:
1581: suffix: 1
1582: args: -no_view
1583: requires: !complex double
1584: filter: grep -v atomic
1585: filter_output: grep -v atomic
1586: test:
1587: suffix: 1_matis
1588: requires: !complex double
1589: args: -no_view -dm_mat_type is
1590: filter: grep -v atomic
1591: filter_output: grep -v atomic
1592: testset:
1593: nsize: 4
1594: requires: !complex double
1595: 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
1596: filter: grep -v atomic
1597: filter_output: grep -v atomic
1598: test:
1599: suffix: fetidp
1600: args: -stokes_fetidp_bddc_pc_bddc_coarse_redundant_pc_type svd
1601: test:
1602: suffix: fetidp_lumped
1603: 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
1604: test:
1605: suffix: fetidp_saddlepoint
1606: 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
1607: test:
1608: suffix: fetidp_saddlepoint_lumped
1609: 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
1610: TEST*/