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";
27
28 /* Contributed by Dave May */
29
30 #include <petscksp.h>
31 #include <petscdm.h>
32 #include <petscdmda.h>
33 #include <petscdmswarm.h>
34 #include <petsc/private/dmimpl.h>
35
36 static PetscErrorCode DMDAApplyBoundaryConditions(DM,Mat,Vec);
37
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
43
44
EvaluateBasis_Q1(PetscScalar _xi[],PetscScalar N[])45 static void EvaluateBasis_Q1(PetscScalar _xi[],PetscScalar N[])
46 {
47 PetscScalar xi = _xi[0];
48 PetscScalar eta = _xi[1];
49
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 }
55
EvaluateBasisDerivatives_Q1(PetscScalar _xi[],PetscScalar dN[][NODES_PER_EL])56 static void EvaluateBasisDerivatives_Q1(PetscScalar _xi[],PetscScalar dN[][NODES_PER_EL])
57 {
58 PetscScalar xi = _xi[0];
59 PetscScalar eta = _xi[1];
60
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);
65
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 }
71
EvaluateDerivatives(PetscScalar dN[][NODES_PER_EL],PetscScalar dNx[][NODES_PER_EL],PetscScalar coords[],PetscScalar * det_J)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;
77
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];
82
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);
89
90 iJ00 = J11/J;
91 iJ01 = -J01/J;
92 iJ10 = -J10/J;
93 iJ11 = J00/J;
94
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 }
99
100 if (det_J) *det_J = J;
101 }
102
CreateGaussQuadrature(PetscInt * ngp,PetscScalar gp_xi[][2],PetscScalar gp_weight[])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 }
115
DMDAGetElementEqnums_up(const PetscInt element[],PetscInt s_u[],PetscInt s_p[])116 static PetscErrorCode DMDAGetElementEqnums_up(const PetscInt element[],PetscInt s_u[],PetscInt s_p[])
117 {
118 PetscInt i;
119 PetscFunctionBeginUser;
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 PetscFunctionReturn(0);
128 }
129
map_wIwDI_uJuDJ(PetscInt wi,PetscInt wd,PetscInt w_NPE,PetscInt w_dof,PetscInt ui,PetscInt ud,PetscInt u_NPE,PetscInt u_dof)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;
133
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 }
140
BForm_DivT(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])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;
148
149 /* define quadrature rule */
150 CreateGaussQuadrature(&ngp,gp_xi,gp_weight);
151
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);
156
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];
160
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 }
165
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];
169
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 }
185
BForm_Grad(PetscScalar Ke[],PetscScalar coords[])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;
192
193 /* define quadrature rule */
194 CreateGaussQuadrature(&ngp,gp_xi,gp_weight);
195
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;
202
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);
208
209 Ke[IJ] -= GNx_p[di][i]*Ni_p[j]*fac;
210 }
211 }
212 }
213 }
214 }
215
BForm_Div(PetscScalar De[],PetscScalar coords[])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;
220
221 PetscMemzero(Ge,sizeof(Ge));
222 BForm_Grad(Ge,coords);
223
224 nr_g = U_DOFS*NODES_PER_EL;
225 nc_g = P_DOFS*NODES_PER_EL;
226
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 }
233
BForm_Stabilisation(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])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;
240
241 /* define quadrature rule */
242 CreateGaussQuadrature(&ngp,gp_xi,gp_weight);
243
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;
250
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 }
257
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 }
269
BForm_ScaledMassMatrix(PetscScalar Ke[],PetscScalar coords[],PetscScalar eta[])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;
276
277 /* define quadrature rule */
278 CreateGaussQuadrature(&ngp,gp_xi,gp_weight);
279
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;
286
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 }
293
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 }
305
LForm_MomentumRHS(PetscScalar Fe[],PetscScalar coords[],PetscScalar fx[],PetscScalar fy[])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;
312
313 /* define quadrature rule */
314 CreateGaussQuadrature(&ngp,gp_xi,gp_weight);
315
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;
322
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 }
329
GetElementCoords(const PetscScalar _coords[],const PetscInt e2n[],PetscScalar el_coords[])330 static PetscErrorCode GetElementCoords(const PetscScalar _coords[],const PetscInt e2n[],PetscScalar el_coords[])
331 {
332 PetscInt i,d;
333 PetscFunctionBeginUser;
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 PetscFunctionReturn(0);
341 }
342
AssembleStokes_A(Mat A,DM stokes_da,DM quadrature)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;
359
360 PetscFunctionBeginUser;
361 ierr = MatZeroEntries(A);CHKERRQ(ierr);
362 /* setup for coords */
363 ierr = DMGetCoordinateDM(stokes_da,&cda);CHKERRQ(ierr);
364 ierr = DMGetCoordinatesLocal(stokes_da,&coords);CHKERRQ(ierr);
365 ierr = VecGetArrayRead(coords,&_coords);CHKERRQ(ierr);
366
367 /* setup for coefficients */
368 ierr = DMSwarmGetField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);CHKERRQ(ierr);
369
370 ierr = DMDAGetElements(stokes_da,&nel,&npe,&element_list);CHKERRQ(ierr);
371 for (eidx = 0; eidx < nel; eidx++) {
372 const PetscInt *element = &element_list[npe*eidx];
373
374 /* get coords for the element */
375 ierr = GetElementCoords(_coords,element,el_coords);CHKERRQ(ierr);
376
377 /* get coefficients for the element */
378 prop_eta = &q_eta[GAUSS_POINTS * eidx];
379
380 /* initialise element stiffness matrix */
381 ierr = PetscMemzero(Ae,sizeof(Ae));CHKERRQ(ierr);
382 ierr = PetscMemzero(Ge,sizeof(Ge));CHKERRQ(ierr);
383 ierr = PetscMemzero(De,sizeof(De));CHKERRQ(ierr);
384 ierr = PetscMemzero(Ce,sizeof(Ce));CHKERRQ(ierr);
385
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);
391
392 /* insert element matrix into global matrix */
393 ierr = DMDAGetElementEqnums_up(element,u_eqn,p_eqn);CHKERRQ(ierr);
394 ierr = MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);CHKERRQ(ierr);
395 ierr = MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);CHKERRQ(ierr);
396 ierr = MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);CHKERRQ(ierr);
397 ierr = MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);CHKERRQ(ierr);
398 }
399 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
400 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
401
402 ierr = DMSwarmRestoreField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);CHKERRQ(ierr);
403 ierr = VecRestoreArrayRead(coords,&_coords);CHKERRQ(ierr);
404 PetscFunctionReturn(0);
405 }
406
AssembleStokes_PC(Mat A,DM stokes_da,DM quadrature)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;
423
424 PetscFunctionBeginUser;
425 ierr = MatZeroEntries(A);CHKERRQ(ierr);
426 /* setup for coords */
427 ierr = DMGetCoordinateDM(stokes_da,&cda);CHKERRQ(ierr);
428 ierr = DMGetCoordinatesLocal(stokes_da,&coords);CHKERRQ(ierr);
429 ierr = VecGetArrayRead(coords,&_coords);CHKERRQ(ierr);
430
431 /* setup for coefficients */
432 ierr = DMSwarmGetField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);CHKERRQ(ierr);
433
434 ierr = DMDAGetElements(stokes_da,&nel,&npe,&element_list);CHKERRQ(ierr);
435 for (eidx = 0; eidx < nel; eidx++) {
436 const PetscInt *element = &element_list[npe*eidx];
437
438 /* get coords for the element */
439 ierr = GetElementCoords(_coords,element,el_coords);CHKERRQ(ierr);
440
441 /* get coefficients for the element */
442 prop_eta = &q_eta[GAUSS_POINTS * eidx];
443
444 /* initialise element stiffness matrix */
445 ierr = PetscMemzero(Ae,sizeof(Ae));CHKERRQ(ierr);
446 ierr = PetscMemzero(Ge,sizeof(Ge));CHKERRQ(ierr);
447 ierr = PetscMemzero(De,sizeof(De));CHKERRQ(ierr);
448 ierr = PetscMemzero(Ce,sizeof(Ce));CHKERRQ(ierr);
449
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);
454
455 /* insert element matrix into global matrix */
456 ierr = DMDAGetElementEqnums_up(element,u_eqn,p_eqn);CHKERRQ(ierr);
457 ierr = MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*U_DOFS,u_eqn,Ae,ADD_VALUES);CHKERRQ(ierr);
458 ierr = MatSetValuesLocal(A,NODES_PER_EL*U_DOFS,u_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ge,ADD_VALUES);CHKERRQ(ierr);
459 ierr = MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*U_DOFS,u_eqn,De,ADD_VALUES);CHKERRQ(ierr);
460 ierr = MatSetValuesLocal(A,NODES_PER_EL*P_DOFS,p_eqn,NODES_PER_EL*P_DOFS,p_eqn,Ce,ADD_VALUES);CHKERRQ(ierr);
461 }
462 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
463 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
464
465 ierr = DMSwarmRestoreField(quadrature,"eta_q",NULL,NULL,(void**)&q_eta);CHKERRQ(ierr);
466 ierr = VecRestoreArrayRead(coords,&_coords);CHKERRQ(ierr);
467
468 PetscFunctionReturn(0);
469 }
470
AssembleStokes_RHS(Vec F,DM stokes_da,DM quadrature)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;
487
488 PetscFunctionBeginUser;
489 ierr = VecZeroEntries(F);CHKERRQ(ierr);
490 /* setup for coords */
491 ierr = DMGetCoordinateDM(stokes_da,&cda);CHKERRQ(ierr);
492 ierr = DMGetCoordinatesLocal(stokes_da,&coords);CHKERRQ(ierr);
493 ierr = VecGetArrayRead(coords,&_coords);CHKERRQ(ierr);
494
495 /* setup for coefficients */
496 ierr = DMSwarmGetField(quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);CHKERRQ(ierr);
497
498 /* get acces to the vector */
499 ierr = DMGetLocalVector(stokes_da,&local_F);CHKERRQ(ierr);
500 ierr = VecZeroEntries(local_F);CHKERRQ(ierr);
501 ierr = VecGetArray(local_F,&LA_F);CHKERRQ(ierr);
502
503 ierr = DMDAGetElements(stokes_da,&nel,&npe,&element_list);CHKERRQ(ierr);
504 for (eidx = 0; eidx < nel; eidx++) {
505 const PetscInt *element = &element_list[npe*eidx];
506
507 /* get coords for the element */
508 ierr = GetElementCoords(_coords,element,el_coords);CHKERRQ(ierr);
509
510 /* get coefficients for the element */
511 prop_fy = &q_rhs[GAUSS_POINTS * eidx];
512
513 /* initialise element stiffness matrix */
514 ierr = PetscMemzero(Fe,sizeof(Fe));CHKERRQ(ierr);
515 ierr = PetscMemzero(He,sizeof(He));CHKERRQ(ierr);
516
517 /* form element stiffness matrix */
518 LForm_MomentumRHS(Fe,el_coords,NULL,prop_fy);
519
520 /* insert element matrix into global matrix */
521 ierr = DMDAGetElementEqnums_up(element,u_eqn,p_eqn);CHKERRQ(ierr);
522
523 for (i=0; i<NODES_PER_EL*U_DOFS; i++) {
524 LA_F[ u_eqn[i] ] += Fe[i];
525 }
526 }
527 ierr = DMSwarmRestoreField(quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);CHKERRQ(ierr);
528 ierr = VecRestoreArrayRead(coords,&_coords);CHKERRQ(ierr);
529
530 ierr = VecRestoreArray(local_F,&LA_F);CHKERRQ(ierr);
531 ierr = DMLocalToGlobalBegin(stokes_da,local_F,ADD_VALUES,F);CHKERRQ(ierr);
532 ierr = DMLocalToGlobalEnd(stokes_da,local_F,ADD_VALUES,F);CHKERRQ(ierr);
533 ierr = DMRestoreLocalVector(stokes_da,&local_F);CHKERRQ(ierr);
534
535 PetscFunctionReturn(0);
536 }
537
DMSwarmPICInsertPointsCellwise(DM dm,DM dmc,PetscInt e,PetscInt npoints,PetscReal xi[],PetscBool proximity_initialization)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;
548
549 PetscFunctionBeginUser;
550 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
551 ierr = DMDAGetElements(dmc,&nel,&npe,&element_list);CHKERRQ(ierr);
552
553 ierr = PetscMalloc1(dim*npoints,&xp);CHKERRQ(ierr);
554 ierr = PetscMalloc1(dim*npe,&elcoor);CHKERRQ(ierr);
555 ierr = PetscMalloc1(npoints,&basis);CHKERRQ(ierr);
556 for (q=0; q<npoints; q++) {
557 ierr = PetscMalloc1(npe,&basis[q]);CHKERRQ(ierr);
558
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;
570
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 }
583
584 ierr = DMGetCoordinatesLocal(dmc,&coor);CHKERRQ(ierr);
585 ierr = VecGetArrayRead(coor,&_coor);CHKERRQ(ierr);
586 /* compute and store the coordinates for the new points */
587 {
588 const PetscInt *element = &element_list[npe*e];
589
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 ierr = VecRestoreArrayRead(coor,&_coor);CHKERRQ(ierr);
607 ierr = DMDARestoreElements(dmc,&nel,&npe,&element_list);CHKERRQ(ierr);
608
609 ierr = DMSwarmGetLocalSize(dm,&ncurr);CHKERRQ(ierr);
610 ierr = DMSwarmAddNPoints(dm,npoints);CHKERRQ(ierr);
611
612 if (proximity_initialization) {
613 PetscInt *nnlist;
614 PetscReal *coor_q,*coor_qn;
615 PetscInt npoints_e,*plist_e;
616
617 ierr = DMSwarmSortGetPointsPerCell(dm,e,&npoints_e,&plist_e);CHKERRQ(ierr);
618
619 ierr = PetscMalloc1(npoints,&nnlist);CHKERRQ(ierr);
620 /* find nearest neighour points in this cell */
621 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
622 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
623 for (q=0; q<npoints; q++) {
624 PetscInt qn,nearest_neighour = -1;
625 PetscReal sep,min_sep = PETSC_MAX_REAL;
626
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 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
643 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
644
645 /* copies the nearest neighbour (nnlist[q]) into the new slot (ncurr+q) */
646 for (q=0; q<npoints; q++) {
647 ierr = DMSwarmCopyPoint(dm,nnlist[q],ncurr+q);CHKERRQ(ierr);
648 }
649 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
650 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
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 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
660 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
661
662 ierr = PetscFree(plist_e);CHKERRQ(ierr);
663 ierr = PetscFree(nnlist);CHKERRQ(ierr);
664 } else {
665 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
666 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
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 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
676 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
677 }
678
679 ierr = PetscFree(xp);CHKERRQ(ierr);
680 ierr = PetscFree(elcoor);CHKERRQ(ierr);
681 for (q=0; q<npoints; q++) {
682 ierr = PetscFree(basis[q]);CHKERRQ(ierr);
683 }
684 ierr = PetscFree(basis);CHKERRQ(ierr);
685 PetscFunctionReturn(0);
686 }
687
MaterialPoint_PopulateCell(DM dm_vp,DM dm_mpoint)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;
697
698 PetscFunctionBeginUser;
699 ierr = DMDAGetElements(dm_vp,&_nel,&_npe,&element);CHKERRQ(ierr);
700 nel = _nel;
701 ierr = DMDARestoreElements(dm_vp,&_nel,&_npe,&element);CHKERRQ(ierr);
702
703 ierr = PetscDTGaussTensorQuadrature(2,1,4,-1.0,1.0,&quadrature);CHKERRQ(ierr);
704 ierr = PetscQuadratureGetData(quadrature,NULL,NULL,&npoints_q,&xi,NULL);CHKERRQ(ierr);
705 ierr = DMSwarmGetCellDM(dm_mpoint,&dmc);CHKERRQ(ierr);
706
707 ierr = DMSwarmSortGetAccess(dm_mpoint);CHKERRQ(ierr);
708
709 cnt = 0;
710 for (e=0; e<nel; e++) {
711 PetscInt npoints_per_cell;
712
713 ierr = DMSwarmSortGetNumberOfPointsPerCell(dm_mpoint,e,&npoints_per_cell);CHKERRQ(ierr);
714
715 if (npoints_per_cell < 12) {
716 ierr = DMSwarmPICInsertPointsCellwise(dm_mpoint,dm_vp,e,npoints_q,(PetscReal*)xi,PETSC_TRUE);CHKERRQ(ierr);
717 cnt++;
718 }
719 }
720 ierr = MPI_Allreduce(&cnt,&cnt_g,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRQ(ierr);
721 if (cnt_g > 0) {
722 ierr = PetscPrintf(PETSC_COMM_WORLD,".... ....pop cont: adjusted %D cells\n",cnt_g);CHKERRQ(ierr);
723 }
724
725 ierr = DMSwarmSortRestoreAccess(dm_mpoint);CHKERRQ(ierr);
726 ierr = PetscQuadratureDestroy(&quadrature);CHKERRQ(ierr);
727 PetscFunctionReturn(0);
728 }
729
MaterialPoint_AdvectRK1(DM dm_vp,Vec vp,PetscReal dt,DM dm_mpoint)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];
743
744 PetscFunctionBeginUser;
745 ierr = DMGetCoordinatesLocal(dm_vp,&coor_l);CHKERRQ(ierr);
746 ierr = VecGetArrayRead(coor_l,&LA_coor);CHKERRQ(ierr);
747
748 ierr = DMGetLocalVector(dm_vp,&vp_l);CHKERRQ(ierr);
749 ierr = DMGlobalToLocalBegin(dm_vp,vp,INSERT_VALUES,vp_l);CHKERRQ(ierr);
750 ierr = DMGlobalToLocalEnd(dm_vp,vp,INSERT_VALUES,vp_l);CHKERRQ(ierr);
751 ierr = VecGetArrayRead(vp_l,&LA_vp);CHKERRQ(ierr);
752
753 ierr = DMDAGetElements(dm_vp,&nel,&npe,&element_list);CHKERRQ(ierr);
754 ierr = DMSwarmGetLocalSize(dm_mpoint,&npoints);CHKERRQ(ierr);
755 ierr = DMSwarmGetField(dm_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr);
756 ierr = DMSwarmGetField(dm_mpoint,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr);
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;
762
763 e = mpfield_cell[p];
764 coor_p = &mpfield_coor[NSD*p];
765 element = &element_list[NODES_PER_EL*e];
766
767 /* compute local coordinates: (xp-x0)/dx = (xip+1)/2 */
768 x0 = &LA_coor[NSD*element[0]];
769 x2 = &LA_coor[NSD*element[2]];
770
771 dx[0] = x2[0] - x0[0];
772 dx[1] = x2[1] - x0[1];
773
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);
780
781 /* evaluate basis functions */
782 EvaluateBasis_Q1(xi_p,Ni);
783
784 /* get cell nodal velocities */
785 for (i=0; i<NODES_PER_EL; i++) {
786 PetscInt nid;
787
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 }
792
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 }
799
800 coor_p[0] += dt * PetscRealPart(vel_p[0]);
801 coor_p[1] += dt * PetscRealPart(vel_p[1]);
802 }
803
804 ierr = DMSwarmRestoreField(dm_mpoint,DMSwarmPICField_cellid,NULL,NULL,(void**)&mpfield_cell);CHKERRQ(ierr);
805 ierr = DMSwarmRestoreField(dm_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&mpfield_coor);CHKERRQ(ierr);
806 ierr = DMDARestoreElements(dm_vp,&nel,&npe,&element_list);CHKERRQ(ierr);
807 ierr = VecRestoreArrayRead(vp_l,&LA_vp);CHKERRQ(ierr);
808 ierr = DMRestoreLocalVector(dm_vp,&vp_l);CHKERRQ(ierr);
809 ierr = VecRestoreArrayRead(coor_l,&LA_coor);CHKERRQ(ierr);
810 PetscFunctionReturn(0);
811 }
812
MaterialPoint_Interpolate(DM dm,Vec eta_v,Vec rho_v,DM dm_quadrature)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;
824 PetscErrorCode ierr;
825
826 PetscFunctionBeginUser;
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 }
832
833 ierr = DMGetLocalVector(dm,&eta_l);CHKERRQ(ierr);
834 ierr = DMGetLocalVector(dm,&rho_l);CHKERRQ(ierr);
835
836 ierr = DMGlobalToLocalBegin(dm,eta_v,INSERT_VALUES,eta_l);CHKERRQ(ierr);
837 ierr = DMGlobalToLocalEnd(dm,eta_v,INSERT_VALUES,eta_l);CHKERRQ(ierr);
838 ierr = DMGlobalToLocalBegin(dm,rho_v,INSERT_VALUES,rho_l);CHKERRQ(ierr);
839 ierr = DMGlobalToLocalEnd(dm,rho_v,INSERT_VALUES,rho_l);CHKERRQ(ierr);
840
841 ierr = VecGetArray(eta_l,&_eta_l);CHKERRQ(ierr);
842 ierr = VecGetArray(rho_l,&_rho_l);CHKERRQ(ierr);
843
844 ierr = DMSwarmGetField(dm_quadrature,"eta_q",NULL,NULL,(void**)&q_eta);CHKERRQ(ierr);
845 ierr = DMSwarmGetField(dm_quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);CHKERRQ(ierr);
846
847 ierr = DMDAGetElements(dm,&nel,&npe,&element_list);CHKERRQ(ierr);
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];
852
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 }
857
858 for (q=0; q<nqp; q++) {
859 PetscScalar eta_q,rho_q;
860
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 }
866
867 q_eta[nqp*e+q] = PetscRealPart(eta_q);
868 q_rhs[nqp*e+q] = PetscRealPart(rho_q);
869 }
870 }
871 ierr = DMDARestoreElements(dm,&nel,&npe,&element_list);CHKERRQ(ierr);
872
873 ierr = DMSwarmRestoreField(dm_quadrature,"rho_q",NULL,NULL,(void**)&q_rhs);CHKERRQ(ierr);
874 ierr = DMSwarmRestoreField(dm_quadrature,"eta_q",NULL,NULL,(void**)&q_eta);CHKERRQ(ierr);
875
876 ierr = VecRestoreArray(rho_l,&_rho_l);CHKERRQ(ierr);
877 ierr = VecRestoreArray(eta_l,&_eta_l);CHKERRQ(ierr);
878 ierr = DMRestoreLocalVector(dm,&rho_l);CHKERRQ(ierr);
879 ierr = DMRestoreLocalVector(dm,&eta_l);CHKERRQ(ierr);
880 PetscFunctionReturn(0);
881 }
882
SolveTimeDepStokes(PetscInt mx,PetscInt my)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;
909
910 PetscFunctionBeginUser;
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 ierr = 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);CHKERRQ(ierr);
922 ierr = DMDASetElementType(dm_stokes,DMDA_ELEMENT_Q1);CHKERRQ(ierr);
923 ierr = DMSetMatType(dm_stokes,MATAIJ);CHKERRQ(ierr);
924 ierr = DMSetFromOptions(dm_stokes);CHKERRQ(ierr);
925 ierr = DMSetUp(dm_stokes);CHKERRQ(ierr);
926 ierr = DMDASetFieldName(dm_stokes,0,"ux");CHKERRQ(ierr);
927 ierr = DMDASetFieldName(dm_stokes,1,"uy");CHKERRQ(ierr);
928 ierr = DMDASetFieldName(dm_stokes,2,"p");CHKERRQ(ierr);
929
930 /* unit box [0,0.9142] x [0,1] */
931 ierr = DMDASetUniformCoordinates(dm_stokes,0.0,0.9142,0.0,1.0,0.,0.);CHKERRQ(ierr);
932 dh = 1.0/((PetscReal)(mx));
933
934 /* Get local number of elements */
935 {
936 ierr = DMDAGetElements(dm_stokes,&nel,&npe,&element_list);CHKERRQ(ierr);
937
938 nel_local = nel;
939
940 ierr = DMDARestoreElements(dm_stokes,&nel,&npe,&element_list);CHKERRQ(ierr);
941 }
942
943 /* Create DMDA for representing scalar fields */
944 ierr = DMDACreateCompatibleDMDA(dm_stokes,1,&dm_coeff);CHKERRQ(ierr);
945
946 /* Create the swarm for storing quadrature point values */
947 ierr = DMCreate(PETSC_COMM_WORLD,&dms_quadrature);CHKERRQ(ierr);
948 ierr = DMSetType(dms_quadrature,DMSWARM);CHKERRQ(ierr);
949 ierr = DMSetDimension(dms_quadrature,2);CHKERRQ(ierr);
950
951 /* Register fields for viscosity and density on the quadrature points */
952 ierr = DMSwarmRegisterPetscDatatypeField(dms_quadrature,"eta_q",1,PETSC_REAL);CHKERRQ(ierr);
953 ierr = DMSwarmRegisterPetscDatatypeField(dms_quadrature,"rho_q",1,PETSC_REAL);CHKERRQ(ierr);
954 ierr = DMSwarmFinalizeFieldRegister(dms_quadrature);CHKERRQ(ierr);
955 ierr = DMSwarmSetLocalSizes(dms_quadrature,nel_local * GAUSS_POINTS,0);CHKERRQ(ierr);
956
957 /* Create the material point swarm */
958 ierr = DMCreate(PETSC_COMM_WORLD,&dms_mpoint);CHKERRQ(ierr);
959 ierr = DMSetType(dms_mpoint,DMSWARM);CHKERRQ(ierr);
960 ierr = DMSetDimension(dms_mpoint,2);CHKERRQ(ierr);
961
962 /* Configure the material point swarm to be of type Particle-In-Cell */
963 ierr = DMSwarmSetType(dms_mpoint,DMSWARM_PIC);CHKERRQ(ierr);
964
965 /*
966 Specify the DM to use for point location and projections
967 within the context of a PIC scheme
968 */
969 ierr = DMSwarmSetCellDM(dms_mpoint,dm_coeff);CHKERRQ(ierr);
970
971 /* Register fields for viscosity and density */
972 ierr = DMSwarmRegisterPetscDatatypeField(dms_mpoint,"eta",1,PETSC_REAL);CHKERRQ(ierr);
973 ierr = DMSwarmRegisterPetscDatatypeField(dms_mpoint,"rho",1,PETSC_REAL);CHKERRQ(ierr);
974 ierr = DMSwarmFinalizeFieldRegister(dms_mpoint);CHKERRQ(ierr);
975
976 ierr = PetscOptionsGetInt(NULL,NULL,"-ppcell",&ppcell,NULL);CHKERRQ(ierr);
977 ierr = DMSwarmSetLocalSizes(dms_mpoint,nel_local * ppcell,100);CHKERRQ(ierr);
978
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 ierr = DMSwarmInsertPointsUsingCellDM(dms_mpoint,DMSWARMPIC_LAYOUT_SUBDIVISION,ppcell);CHKERRQ(ierr);
993
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];
1000
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 ierr = DMSwarmSetPointsUniformCoordinates(dms_mpoint,min,max,npoints_dir_x,ADD_VALUES);CHKERRQ(ierr);
1006 }
1007
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];
1016
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 ierr = DMSwarmSetPointsUniformCoordinates(dms_mpoint,min,max,npoints_dir_x,ADD_VALUES);CHKERRQ(ierr);
1022 }
1023
1024 ierr = DMView(dms_mpoint,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
1025
1026 /* Define initial material properties on each particle in the material point swarm */
1027 ierr = PetscOptionsGetReal(NULL,NULL,"-delta_eta",&delta_eta,NULL);CHKERRQ(ierr);
1028 ierr = PetscOptionsGetBool(NULL,NULL,"-randomize_coords",&randomize_coords,NULL);CHKERRQ(ierr);
1029 ierr = PetscOptionsGetReal(NULL,NULL,"-randomize_fac",&randomize_fac,NULL);CHKERRQ(ierr);
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;
1036
1037 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr);
1038
1039 ierr = PetscRandomCreate(PETSC_COMM_SELF,&r);CHKERRQ(ierr);
1040 ierr = PetscRandomSetInterval(r,-randomize_fac*dh,randomize_fac*dh);CHKERRQ(ierr);
1041 ierr = PetscRandomSetSeed(r,(unsigned long)rank);CHKERRQ(ierr);
1042 ierr = PetscRandomSeed(r);CHKERRQ(ierr);
1043
1044 ierr = DMDAGetElements(dm_stokes,&nel,&npe,&element_list);CHKERRQ(ierr);
1045
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 ierr = DMSwarmGetField(dms_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
1056 ierr = DMSwarmGetField(dms_mpoint,"eta", NULL,NULL,(void**)&array_e);CHKERRQ(ierr);
1057 ierr = DMSwarmGetField(dms_mpoint,"rho", NULL,NULL,(void**)&array_r);CHKERRQ(ierr);
1058
1059 ierr = DMSwarmGetLocalSize(dms_mpoint,&npoints);CHKERRQ(ierr);
1060 for (p = 0; p < npoints; p++) {
1061 PetscReal x_p[2],rr[2];
1062
1063 if (randomize_coords) {
1064 ierr = PetscRandomGetValueReal(r,&rr[0]);CHKERRQ(ierr);
1065 ierr = PetscRandomGetValueReal(r,&rr[1]);CHKERRQ(ierr);
1066 array_x[2*p + 0] += rr[0];
1067 array_x[2*p + 1] += rr[1];
1068 }
1069
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];
1073
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 }
1084
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 ierr = DMSwarmRestoreField(dms_mpoint,"rho",NULL,NULL,(void**)&array_r);CHKERRQ(ierr);
1091 ierr = DMSwarmRestoreField(dms_mpoint,"eta",NULL,NULL,(void**)&array_e);CHKERRQ(ierr);
1092 ierr = DMSwarmRestoreField(dms_mpoint,DMSwarmPICField_coor,NULL,NULL,(void**)&array_x);CHKERRQ(ierr);
1093
1094 ierr = DMDARestoreElements(dm_stokes,&nel,&npe,&element_list);CHKERRQ(ierr);
1095 ierr = PetscRandomDestroy(&r);CHKERRQ(ierr);
1096 }
1097
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 ierr = DMSwarmMigrate(dms_mpoint,PETSC_TRUE);CHKERRQ(ierr);
1104 }
1105
1106 ierr = PetscOptionsGetBool(NULL,NULL,"-no_view",&no_view,NULL);CHKERRQ(ierr);
1107 if (!no_view) {
1108 ierr = DMSwarmViewXDMF(dms_mpoint,"ic_coeff_dms.xmf");CHKERRQ(ierr);
1109 }
1110
1111 /* project the swarm properties */
1112 ierr = DMSwarmProjectFields(dms_mpoint,2,fieldnames,&pfields,PETSC_FALSE);CHKERRQ(ierr);
1113 eta_v = pfields[0];
1114 rho_v = pfields[1];
1115 ierr = PetscObjectSetName((PetscObject)eta_v,"eta");CHKERRQ(ierr);
1116 ierr = PetscObjectSetName((PetscObject)rho_v,"rho");CHKERRQ(ierr);
1117 ierr = MaterialPoint_Interpolate(dm_coeff,eta_v,rho_v,dms_quadrature);CHKERRQ(ierr);
1118
1119 /* view projected coefficients eta and rho */
1120 if (!no_view) {
1121 PetscViewer viewer;
1122
1123 ierr = PetscViewerCreate(PETSC_COMM_WORLD,&viewer);CHKERRQ(ierr);
1124 ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr);
1125 ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
1126 ierr = PetscViewerFileSetName(viewer,"ic_coeff_dmda.vts");CHKERRQ(ierr);
1127 ierr = VecView(eta_v,viewer);CHKERRQ(ierr);
1128 ierr = VecView(rho_v,viewer);CHKERRQ(ierr);
1129 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
1130 }
1131
1132 ierr = DMCreateMatrix(dm_stokes,&A);CHKERRQ(ierr);
1133 ierr = DMCreateMatrix(dm_stokes,&B);CHKERRQ(ierr);
1134 ierr = DMCreateGlobalVector(dm_stokes,&f);CHKERRQ(ierr);
1135 ierr = DMCreateGlobalVector(dm_stokes,&X);CHKERRQ(ierr);
1136
1137 ierr = AssembleStokes_A(A,dm_stokes,dms_quadrature);CHKERRQ(ierr);
1138 ierr = AssembleStokes_PC(B,dm_stokes,dms_quadrature);CHKERRQ(ierr);
1139 ierr = AssembleStokes_RHS(f,dm_stokes,dms_quadrature);CHKERRQ(ierr);
1140
1141 ierr = DMDAApplyBoundaryConditions(dm_stokes,A,f);CHKERRQ(ierr);
1142 ierr = DMDAApplyBoundaryConditions(dm_stokes,B,NULL);CHKERRQ(ierr);
1143
1144 ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
1145 ierr = KSPSetOptionsPrefix(ksp,"stokes_");CHKERRQ(ierr);
1146 ierr = KSPSetDM(ksp,dm_stokes);CHKERRQ(ierr);
1147 ierr = KSPSetDMActive(ksp,PETSC_FALSE);CHKERRQ(ierr);
1148 ierr = KSPSetOperators(ksp,A,B);CHKERRQ(ierr);
1149 ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
1150 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
1151 ierr = PetscObjectTypeCompare((PetscObject)pc,PCBDDC,&isbddc);CHKERRQ(ierr);
1152 if (isbddc) {
1153 ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
1154 }
1155
1156 /* Define u-v-p indices for fieldsplit */
1157 {
1158 PC pc;
1159 const PetscInt ufields[] = {0,1},pfields[1] = {2};
1160
1161 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
1162 ierr = PCFieldSplitSetBlockSize(pc,3);CHKERRQ(ierr);
1163 ierr = PCFieldSplitSetFields(pc,"u",2,ufields,ufields);CHKERRQ(ierr);
1164 ierr = PCFieldSplitSetFields(pc,"p",1,pfields,pfields);CHKERRQ(ierr);
1165 }
1166
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;
1174
1175 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
1176
1177 is_pcfs = PETSC_FALSE;
1178 ierr = PetscObjectTypeCompare((PetscObject)pc,PCFIELDSPLIT,&is_pcfs);CHKERRQ(ierr);
1179
1180 if (is_pcfs) {
1181 ierr = KSPSetUp(ksp);CHKERRQ(ierr);
1182 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
1183 ierr = PCFieldSplitGetSubKSP(pc,&nsplits,&sub_ksp);CHKERRQ(ierr);
1184 ksp_u = sub_ksp[0];
1185 ierr = PetscFree(sub_ksp);CHKERRQ(ierr);
1186
1187 if (nsplits == 2) {
1188 ierr = DMDACreateCompatibleDMDA(dm_stokes,2,&dm_u);CHKERRQ(ierr);
1189
1190 ierr = KSPSetDM(ksp_u,dm_u);CHKERRQ(ierr);
1191 ierr = KSPSetDMActive(ksp_u,PETSC_FALSE);CHKERRQ(ierr);
1192 ierr = DMDestroy(&dm_u);CHKERRQ(ierr);
1193
1194 /* enforce galerkin coarse grids be used */
1195 ierr = KSPGetPC(ksp_u,&pc_u);CHKERRQ(ierr);
1196 ierr = PCMGSetGalerkin(pc_u,PC_MG_GALERKIN_PMAT);CHKERRQ(ierr);
1197 }
1198 }
1199 }
1200
1201 dump_freq = 10;
1202 ierr = PetscOptionsGetInt(NULL,NULL,"-dump_freq",&dump_freq,NULL);CHKERRQ(ierr);
1203 nt = 10;
1204 ierr = PetscOptionsGetInt(NULL,NULL,"-nt",&nt,NULL);CHKERRQ(ierr);
1205 time = 0.0;
1206 for (tk=1; tk<=nt; tk++) {
1207
1208 ierr = PetscPrintf(PETSC_COMM_WORLD,".... assemble\n");CHKERRQ(ierr);
1209 ierr = AssembleStokes_A(A,dm_stokes,dms_quadrature);CHKERRQ(ierr);
1210 ierr = AssembleStokes_PC(B,dm_stokes,dms_quadrature);CHKERRQ(ierr);
1211 ierr = AssembleStokes_RHS(f,dm_stokes,dms_quadrature);CHKERRQ(ierr);
1212
1213 ierr = PetscPrintf(PETSC_COMM_WORLD,".... bc imposition\n");CHKERRQ(ierr);
1214 ierr = DMDAApplyBoundaryConditions(dm_stokes,A,f);CHKERRQ(ierr);
1215 ierr = DMDAApplyBoundaryConditions(dm_stokes,B,NULL);CHKERRQ(ierr);
1216
1217 ierr = PetscPrintf(PETSC_COMM_WORLD,".... solve\n");CHKERRQ(ierr);
1218 ierr = KSPSetOperators(ksp,A, isbddc ? A : B);CHKERRQ(ierr);
1219 ierr = KSPSolve(ksp,f,X);CHKERRQ(ierr);
1220
1221 ierr = VecStrideMax(X,0,NULL,&vx[1]);CHKERRQ(ierr);
1222 ierr = VecStrideMax(X,1,NULL,&vy[1]);CHKERRQ(ierr);
1223 ierr = VecStrideMin(X,0,NULL,&vx[0]);CHKERRQ(ierr);
1224 ierr = VecStrideMin(X,1,NULL,&vy[0]);CHKERRQ(ierr);
1225
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);
1230
1231 dt_max = 2.0;
1232 dt = 0.5 * (dh / max_v_step);
1233 ierr = 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);CHKERRQ(ierr);
1234 dt = PetscMin(dt_max,dt);
1235
1236 /* advect */
1237 ierr = PetscPrintf(PETSC_COMM_WORLD,".... advect\n");CHKERRQ(ierr);
1238 ierr = MaterialPoint_AdvectRK1(dm_stokes,X,dt,dms_mpoint);CHKERRQ(ierr);
1239
1240 /* migrate */
1241 ierr = PetscPrintf(PETSC_COMM_WORLD,".... migrate\n");CHKERRQ(ierr);
1242 ierr = DMSwarmMigrate(dms_mpoint,PETSC_TRUE);CHKERRQ(ierr);
1243
1244 /* update cell population */
1245 ierr = PetscPrintf(PETSC_COMM_WORLD,".... populate cells\n");CHKERRQ(ierr);
1246 ierr = MaterialPoint_PopulateCell(dm_stokes,dms_mpoint);CHKERRQ(ierr);
1247
1248 /* update coefficients on quadrature points */
1249 ierr = PetscPrintf(PETSC_COMM_WORLD,".... project\n");CHKERRQ(ierr);
1250 ierr = DMSwarmProjectFields(dms_mpoint,2,fieldnames,&pfields,PETSC_TRUE);CHKERRQ(ierr);
1251 eta_v = pfields[0];
1252 rho_v = pfields[1];
1253 ierr = PetscPrintf(PETSC_COMM_WORLD,".... interp\n");CHKERRQ(ierr);
1254 ierr = MaterialPoint_Interpolate(dm_coeff,eta_v,rho_v,dms_quadrature);CHKERRQ(ierr);
1255
1256 if (tk%dump_freq == 0) {
1257 PetscViewer viewer;
1258
1259 ierr = PetscPrintf(PETSC_COMM_WORLD,".... write XDMF, VTS\n");CHKERRQ(ierr);
1260 ierr = PetscSNPrintf(filename,PETSC_MAX_PATH_LEN-1,"step%.4D_coeff_dms.xmf",tk);CHKERRQ(ierr);
1261 ierr = DMSwarmViewXDMF(dms_mpoint,filename);CHKERRQ(ierr);
1262
1263 ierr = PetscSNPrintf(filename,PETSC_MAX_PATH_LEN-1,"step%.4D_vp_dm.vts",tk);CHKERRQ(ierr);
1264 ierr = PetscViewerCreate(PETSC_COMM_WORLD,&viewer);CHKERRQ(ierr);
1265 ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr);
1266 ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
1267 ierr = PetscViewerFileSetName(viewer,filename);CHKERRQ(ierr);
1268 ierr = VecView(X,viewer);CHKERRQ(ierr);
1269 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
1270 }
1271 time += dt;
1272 ierr = PetscPrintf(PETSC_COMM_WORLD,"step %D : time %1.2e \n",tk,(double)time);CHKERRQ(ierr);
1273 }
1274
1275 ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
1276 ierr = VecDestroy(&X);CHKERRQ(ierr);
1277 ierr = VecDestroy(&f);CHKERRQ(ierr);
1278 ierr = MatDestroy(&A);CHKERRQ(ierr);
1279 ierr = MatDestroy(&B);CHKERRQ(ierr);
1280 ierr = VecDestroy(&eta_v);CHKERRQ(ierr);
1281 ierr = VecDestroy(&rho_v);CHKERRQ(ierr);
1282 ierr = PetscFree(pfields);CHKERRQ(ierr);
1283
1284 ierr = DMDestroy(&dms_mpoint);CHKERRQ(ierr);
1285 ierr = DMDestroy(&dms_quadrature);CHKERRQ(ierr);
1286 ierr = DMDestroy(&dm_coeff);CHKERRQ(ierr);
1287 ierr = DMDestroy(&dm_stokes);CHKERRQ(ierr);
1288 PetscFunctionReturn(0);
1289 }
1290
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 */
main(int argc,char ** args)1295 int main(int argc,char **args)
1296 {
1297 PetscErrorCode ierr;
1298 PetscInt mx,my;
1299 PetscBool set = PETSC_FALSE;
1300
1301 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
1302 mx = my = 10;
1303 ierr = PetscOptionsGetInt(NULL,NULL,"-mx",&mx,NULL);CHKERRQ(ierr);
1304 ierr = PetscOptionsGetInt(NULL,NULL,"-my",&my,NULL);CHKERRQ(ierr);
1305 ierr = PetscOptionsGetInt(NULL,NULL,"-mxy",&mx,&set);CHKERRQ(ierr);
1306 if (set) {
1307 my = mx;
1308 }
1309 ierr = SolveTimeDepStokes(mx,my);CHKERRQ(ierr);
1310 ierr = PetscFinalize();
1311 return ierr;
1312 }
1313
1314 /* -------------------------- helpers for boundary conditions -------------------------------- */
BCApplyZero_EAST(DM da,PetscInt d_idx,Mat A,Vec b)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;
1329
1330 PetscFunctionBeginUser;
1331 ierr = DMGetLocalToGlobalMapping(da,<ogm);CHKERRQ(ierr);
1332 ierr = ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);CHKERRQ(ierr);
1333
1334 ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
1335 ierr = DMGetCoordinatesLocal(da,&coords);CHKERRQ(ierr);
1336 ierr = DMDAVecGetArray(cda,coords,&_coords);CHKERRQ(ierr);
1337 ierr = DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);CHKERRQ(ierr);
1338 ierr = DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);CHKERRQ(ierr);
1339
1340 ierr = PetscMalloc1(ny*n_dofs,&bc_global_ids);CHKERRQ(ierr);
1341 ierr = PetscMalloc1(ny*n_dofs,&bc_vals);CHKERRQ(ierr);
1342
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;
1345
1346 i = nx-1;
1347 for (j = 0; j < ny; j++) {
1348 PetscInt local_id;
1349
1350 local_id = i+j*nx;
1351
1352 bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1353
1354 bc_vals[j] = 0.0;
1355 }
1356 ierr = ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);CHKERRQ(ierr);
1357 nbcs = 0;
1358 if ((si+nx) == (M)) nbcs = ny;
1359
1360 if (b) {
1361 ierr = VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);CHKERRQ(ierr);
1362 ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
1363 ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
1364 }
1365 if (A) {
1366 ierr = MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);CHKERRQ(ierr);
1367 }
1368
1369 ierr = PetscFree(bc_vals);CHKERRQ(ierr);
1370 ierr = PetscFree(bc_global_ids);CHKERRQ(ierr);
1371
1372 ierr = DMDAVecRestoreArray(cda,coords,&_coords);CHKERRQ(ierr);
1373 PetscFunctionReturn(0);
1374 }
1375
BCApplyZero_WEST(DM da,PetscInt d_idx,Mat A,Vec b)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;
1390
1391 PetscFunctionBeginUser;
1392 ierr = DMGetLocalToGlobalMapping(da,<ogm);CHKERRQ(ierr);
1393 ierr = ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);CHKERRQ(ierr);
1394
1395 ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
1396 ierr = DMGetCoordinatesLocal(da,&coords);CHKERRQ(ierr);
1397 ierr = DMDAVecGetArray(cda,coords,&_coords);CHKERRQ(ierr);
1398 ierr = DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);CHKERRQ(ierr);
1399 ierr = DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);CHKERRQ(ierr);
1400
1401 ierr = PetscMalloc1(ny*n_dofs,&bc_global_ids);CHKERRQ(ierr);
1402 ierr = PetscMalloc1(ny*n_dofs,&bc_vals);CHKERRQ(ierr);
1403
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;
1406
1407 i = 0;
1408 for (j = 0; j < ny; j++) {
1409 PetscInt local_id;
1410
1411 local_id = i+j*nx;
1412
1413 bc_global_ids[j] = g_idx[n_dofs*local_id+d_idx];
1414
1415 bc_vals[j] = 0.0;
1416 }
1417 ierr = ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);CHKERRQ(ierr);
1418 nbcs = 0;
1419 if (si == 0) nbcs = ny;
1420
1421 if (b) {
1422 ierr = VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);CHKERRQ(ierr);
1423 ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
1424 ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
1425 }
1426
1427 if (A) {
1428 ierr = MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);CHKERRQ(ierr);
1429 }
1430
1431 ierr = PetscFree(bc_vals);CHKERRQ(ierr);
1432 ierr = PetscFree(bc_global_ids);CHKERRQ(ierr);
1433
1434 ierr = DMDAVecRestoreArray(cda,coords,&_coords);CHKERRQ(ierr);
1435 PetscFunctionReturn(0);
1436 }
1437
BCApplyZero_NORTH(DM da,PetscInt d_idx,Mat A,Vec b)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;
1452
1453 PetscFunctionBeginUser;
1454 ierr = DMGetLocalToGlobalMapping(da,<ogm);CHKERRQ(ierr);
1455 ierr = ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);CHKERRQ(ierr);
1456
1457 ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
1458 ierr = DMGetCoordinatesLocal(da,&coords);CHKERRQ(ierr);
1459 ierr = DMDAVecGetArray(cda,coords,&_coords);CHKERRQ(ierr);
1460 ierr = DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);CHKERRQ(ierr);
1461 ierr = DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);CHKERRQ(ierr);
1462
1463 ierr = PetscMalloc1(nx,&bc_global_ids);CHKERRQ(ierr);
1464 ierr = PetscMalloc1(nx,&bc_vals);CHKERRQ(ierr);
1465
1466 /* init the entries to -1 so VecSetValues will ignore them */
1467 for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1468
1469 j = ny-1;
1470 for (i = 0; i < nx; i++) {
1471 PetscInt local_id;
1472
1473 local_id = i+j*nx;
1474
1475 bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1476
1477 bc_vals[i] = 0.0;
1478 }
1479 ierr = ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);CHKERRQ(ierr);
1480 nbcs = 0;
1481 if ((sj+ny) == (N)) nbcs = nx;
1482
1483 if (b) {
1484 ierr = VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);CHKERRQ(ierr);
1485 ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
1486 ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
1487 }
1488 if (A) {
1489 ierr = MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);CHKERRQ(ierr);
1490 }
1491
1492 ierr = PetscFree(bc_vals);CHKERRQ(ierr);
1493 ierr = PetscFree(bc_global_ids);CHKERRQ(ierr);
1494
1495 ierr = DMDAVecRestoreArray(cda,coords,&_coords);CHKERRQ(ierr);
1496 PetscFunctionReturn(0);
1497 }
1498
BCApplyZero_SOUTH(DM da,PetscInt d_idx,Mat A,Vec b)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;
1513
1514 PetscFunctionBeginUser;
1515 ierr = DMGetLocalToGlobalMapping(da,<ogm);CHKERRQ(ierr);
1516 ierr = ISLocalToGlobalMappingGetIndices(ltogm,&g_idx);CHKERRQ(ierr);
1517
1518 ierr = DMGetCoordinateDM(da,&cda);CHKERRQ(ierr);
1519 ierr = DMGetCoordinatesLocal(da,&coords);CHKERRQ(ierr);
1520 ierr = DMDAVecGetArray(cda,coords,&_coords);CHKERRQ(ierr);
1521 ierr = DMDAGetGhostCorners(cda,&si,&sj,0,&nx,&ny,0);CHKERRQ(ierr);
1522 ierr = DMDAGetInfo(da,0,&M,&N,0,0,0,0,&n_dofs,0,0,0,0,0);CHKERRQ(ierr);
1523
1524 ierr = PetscMalloc1(nx,&bc_global_ids);CHKERRQ(ierr);
1525 ierr = PetscMalloc1(nx,&bc_vals);CHKERRQ(ierr);
1526
1527 /* init the entries to -1 so VecSetValues will ignore them */
1528 for (i = 0; i < nx; i++) bc_global_ids[i] = -1;
1529
1530 j = 0;
1531 for (i = 0; i < nx; i++) {
1532 PetscInt local_id;
1533
1534 local_id = i+j*nx;
1535
1536 bc_global_ids[i] = g_idx[n_dofs*local_id+d_idx];
1537
1538 bc_vals[i] = 0.0;
1539 }
1540 ierr = ISLocalToGlobalMappingRestoreIndices(ltogm,&g_idx);CHKERRQ(ierr);
1541 nbcs = 0;
1542 if (sj == 0) nbcs = nx;
1543
1544 if (b) {
1545 ierr = VecSetValues(b,nbcs,bc_global_ids,bc_vals,INSERT_VALUES);CHKERRQ(ierr);
1546 ierr = VecAssemblyBegin(b);CHKERRQ(ierr);
1547 ierr = VecAssemblyEnd(b);CHKERRQ(ierr);
1548 }
1549 if (A) {
1550 ierr = MatZeroRowsColumns(A,nbcs,bc_global_ids,1.0,0,0);CHKERRQ(ierr);
1551 }
1552
1553 ierr = PetscFree(bc_vals);CHKERRQ(ierr);
1554 ierr = PetscFree(bc_global_ids);CHKERRQ(ierr);
1555
1556 ierr = DMDAVecRestoreArray(cda,coords,&_coords);CHKERRQ(ierr);
1557 PetscFunctionReturn(0);
1558 }
1559
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 */
DMDAApplyBoundaryConditions(DM dm_stokes,Mat A,Vec f)1564 static PetscErrorCode DMDAApplyBoundaryConditions(DM dm_stokes,Mat A,Vec f)
1565 {
1566 PetscErrorCode ierr;
1567
1568 PetscFunctionBeginUser;
1569 ierr = BCApplyZero_NORTH(dm_stokes,0,A,f);CHKERRQ(ierr);
1570 ierr = BCApplyZero_NORTH(dm_stokes,1,A,f);CHKERRQ(ierr);
1571 ierr = BCApplyZero_EAST(dm_stokes,0,A,f);CHKERRQ(ierr);
1572 ierr = BCApplyZero_SOUTH(dm_stokes,0,A,f);CHKERRQ(ierr);
1573 ierr = BCApplyZero_SOUTH(dm_stokes,1,A,f);CHKERRQ(ierr);
1574 ierr = BCApplyZero_WEST(dm_stokes,0,A,f);CHKERRQ(ierr);
1575 PetscFunctionReturn(0);
1576 }
1577
1578 /*TEST
1579
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*/
1611