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