Actual source code: ex63.c
1: static char help[] = "Stokes Problem in 2d and 3d with particles.\n\
2: We solve the Stokes problem in a rectangular\n\
3: domain, using a parallel unstructured mesh (DMPLEX) to discretize it and particles (DMSWARM).\n\n\n";
5: /*
6: The isoviscous Stokes problem, which we discretize using the finite
7: element method on an unstructured mesh. The weak form equations are
9: < \nabla v, \nabla u + {\nabla u}^T > - < \nabla\cdot v, p > + < v, f > = 0
10: < q, \nabla\cdot u > = 0
12: We start with homogeneous Dirichlet conditions.
13: */
15: #include <petscdmplex.h>
16: #include <petscdmswarm.h>
17: #include <petscsnes.h>
18: #include <petscds.h>
20: typedef enum {
21: NEUMANN,
22: DIRICHLET
23: } BCType;
24: typedef enum {
25: RUN_FULL,
26: RUN_TEST
27: } RunType;
29: typedef struct {
30: RunType runType; /* Whether to run tests, or solve the full problem */
31: PetscBool showInitial, showSolution, showError;
32: BCType bcType;
33: } AppCtx;
35: PetscErrorCode zero_scalar(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
36: {
37: u[0] = 0.0;
38: return 0;
39: }
40: PetscErrorCode zero_vector(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
41: {
42: PetscInt d;
43: for (d = 0; d < dim; ++d) u[d] = 0.0;
44: return 0;
45: }
47: PetscErrorCode linear_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
48: {
49: u[0] = x[0];
50: u[1] = 0.0;
51: return 0;
52: }
54: /*
55: In 2D we use exact solution:
57: u = x^2 + y^2
58: v = 2 x^2 - 2xy
59: p = x + y - 1
60: f_x = f_y = 3
62: so that
64: -\Delta u + \nabla p + f = <-4, -4> + <1, 1> + <3, 3> = 0
65: \nabla \cdot u = 2x - 2x = 0
66: */
67: PetscErrorCode quadratic_u_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
68: {
69: u[0] = x[0] * x[0] + x[1] * x[1];
70: u[1] = 2.0 * x[0] * x[0] - 2.0 * x[0] * x[1];
71: return 0;
72: }
74: PetscErrorCode linear_p_2d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx)
75: {
76: *p = x[0] + x[1] - 1.0;
77: return 0;
78: }
79: PetscErrorCode constant_p(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx)
80: {
81: *p = 1.0;
82: return 0;
83: }
85: void f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
86: {
87: PetscInt c;
88: for (c = 0; c < dim; ++c) f0[c] = 3.0;
89: }
91: /* gradU[comp*dim+d] = {u_x, u_y, v_x, v_y} or {u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z}
92: u[Ncomp] = {p} */
93: void f1_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
94: {
95: const PetscInt Ncomp = dim;
96: PetscInt comp, d;
98: for (comp = 0; comp < Ncomp; ++comp) {
99: for (d = 0; d < dim; ++d) {
100: /* f1[comp*dim+d] = 0.5*(gradU[comp*dim+d] + gradU[d*dim+comp]); */
101: f1[comp * dim + d] = u_x[comp * dim + d];
102: }
103: f1[comp * dim + comp] -= u[Ncomp];
104: }
105: }
107: /* gradU[comp*dim+d] = {u_x, u_y, v_x, v_y} or {u_x, u_y, u_z, v_x, v_y, v_z, w_x, w_y, w_z} */
108: void f0_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
109: {
110: PetscInt d;
111: for (d = 0, f0[0] = 0.0; d < dim; ++d) f0[0] += u_x[d * dim + d];
112: }
114: void f1_p(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
115: {
116: PetscInt d;
117: for (d = 0; d < dim; ++d) f1[d] = 0.0;
118: }
120: /* < q, \nabla\cdot u >
121: NcompI = 1, NcompJ = dim */
122: void g1_pu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[])
123: {
124: PetscInt d;
125: for (d = 0; d < dim; ++d) g1[d * dim + d] = 1.0; /* \frac{\partial\phi^{u_d}}{\partial x_d} */
126: }
128: /* -< \nabla\cdot v, p >
129: NcompI = dim, NcompJ = 1 */
130: void g2_up(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[])
131: {
132: PetscInt d;
133: for (d = 0; d < dim; ++d) g2[d * dim + d] = -1.0; /* \frac{\partial\psi^{u_d}}{\partial x_d} */
134: }
136: /* < \nabla v, \nabla u + {\nabla u}^T >
137: This just gives \nabla u, give the perdiagonal for the transpose */
138: void g3_uu(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
139: {
140: const PetscInt Ncomp = dim;
141: PetscInt compI, d;
143: for (compI = 0; compI < Ncomp; ++compI) {
144: for (d = 0; d < dim; ++d) g3[((compI * Ncomp + compI) * dim + d) * dim + d] = 1.0;
145: }
146: }
148: /*
149: In 3D we use exact solution:
151: u = x^2 + y^2
152: v = y^2 + z^2
153: w = x^2 + y^2 - 2(x+y)z
154: p = x + y + z - 3/2
155: f_x = f_y = f_z = 3
157: so that
159: -\Delta u + \nabla p + f = <-4, -4, -4> + <1, 1, 1> + <3, 3, 3> = 0
160: \nabla \cdot u = 2x + 2y - 2(x + y) = 0
161: */
162: PetscErrorCode quadratic_u_3d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
163: {
164: u[0] = x[0] * x[0] + x[1] * x[1];
165: u[1] = x[1] * x[1] + x[2] * x[2];
166: u[2] = x[0] * x[0] + x[1] * x[1] - 2.0 * (x[0] + x[1]) * x[2];
167: return 0;
168: }
170: PetscErrorCode linear_p_3d(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *p, void *ctx)
171: {
172: *p = x[0] + x[1] + x[2] - 1.5;
173: return 0;
174: }
176: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
177: {
178: const char *bcTypes[2] = {"neumann", "dirichlet"};
179: const char *runTypes[2] = {"full", "test"};
180: PetscInt bc, run;
183: options->runType = RUN_FULL;
184: options->bcType = DIRICHLET;
185: options->showInitial = PETSC_FALSE;
186: options->showSolution = PETSC_TRUE;
187: options->showError = PETSC_FALSE;
189: PetscOptionsBegin(comm, "", "Stokes Problem Options", "DMPLEX");
190: run = options->runType;
191: PetscOptionsEList("-run_type", "The run type", "ex62.c", runTypes, 2, runTypes[options->runType], &run, NULL);
192: options->runType = (RunType)run;
193: bc = options->bcType;
194: PetscOptionsEList("-bc_type", "Type of boundary condition", "ex62.c", bcTypes, 2, bcTypes[options->bcType], &bc, NULL);
195: options->bcType = (BCType)bc;
197: PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex62.c", options->showInitial, &options->showInitial, NULL);
198: PetscOptionsBool("-show_solution", "Output the solution for verification", "ex62.c", options->showSolution, &options->showSolution, NULL);
199: PetscOptionsBool("-show_error", "Output the error for verification", "ex62.c", options->showError, &options->showError, NULL);
200: PetscOptionsEnd();
201: return 0;
202: }
204: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
205: {
206: Vec lv;
207: PetscInt p;
208: PetscMPIInt rank, size;
211: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
212: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
213: DMGetLocalVector(dm, &lv);
214: DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv);
215: DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv);
216: PetscPrintf(PETSC_COMM_WORLD, "Local function\n");
217: for (p = 0; p < size; ++p) {
218: if (p == rank) VecView(lv, PETSC_VIEWER_STDOUT_SELF);
219: PetscBarrier((PetscObject)dm);
220: }
221: DMRestoreLocalVector(dm, &lv);
222: return 0;
223: }
225: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
226: {
228: DMCreate(comm, dm);
229: DMSetType(*dm, DMPLEX);
230: DMSetFromOptions(*dm);
231: DMViewFromOptions(*dm, NULL, "-dm_view");
232: return 0;
233: }
235: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
236: {
237: PetscDS prob;
240: DMGetDS(dm, &prob);
241: PetscDSSetResidual(prob, 0, f0_u, f1_u);
242: PetscDSSetResidual(prob, 1, f0_p, f1_p);
243: PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);
244: PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_up, NULL);
245: PetscDSSetJacobian(prob, 1, 0, NULL, g1_pu, NULL, NULL);
246: switch (user->dim) {
247: case 2:
248: user->exactFuncs[0] = quadratic_u_2d;
249: user->exactFuncs[1] = linear_p_2d;
250: break;
251: case 3:
252: user->exactFuncs[0] = quadratic_u_3d;
253: user->exactFuncs[1] = linear_p_3d;
254: break;
255: default:
256: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
257: }
258: return 0;
259: }
261: PetscErrorCode SetupDiscretization(DM dm, AppCtx *user)
262: {
263: DM cdm = dm;
264: const PetscInt dim = user->dim;
265: const PetscInt id = 1;
266: PetscFE fe[2];
267: PetscDS ds;
268: DMLabel label;
269: MPI_Comm comm;
272: /* Create finite element */
273: PetscObjectGetComm((PetscObject)dm, &comm);
274: PetscFECreateDefault(comm, dim, dim, user->simplex, "vel_", PETSC_DEFAULT, &fe[0]);
275: PetscObjectSetName((PetscObject)fe[0], "velocity");
276: PetscFECreateDefault(comm, dim, 1, user->simplex, "pres_", PETSC_DEFAULT, &fe[1]);
277: PetscFECopyQuadrature(fe[0], fe[1]);
278: PetscObjectSetName((PetscObject)fe[1], "pressure");
279: /* Set discretization and boundary conditions for each mesh */
280: while (cdm) {
281: DMGetDS(cdm, &ds);
282: PetscDSSetDiscretization(ds, 0, (PetscObject)fe[0]);
283: PetscDSSetDiscretization(ds, 1, (PetscObject)fe[1]);
284: SetupProblem(cdm, user);
285: if (user->bcType == NEUMANN) DMGetLabel(cdm, "boundary", &label);
286: else DMGetLabel(cdm, "marker", &label);
287: DMAddBoundary(cdm, user->bcType == DIRICHLET ? DM_BC_ESSENTIAL : DM_BC_NATURAL, "wall", label, 1, &id, 0, 0, NULL, (void (*)())user->exactFuncs[0], NULL, user, NULL);
288: DMGetCoarseDM(cdm, &cdm);
289: }
290: PetscFEDestroy(&fe[0]);
291: PetscFEDestroy(&fe[1]);
292: return 0;
293: }
295: PetscErrorCode CreatePressureNullSpace(DM dm, AppCtx *user, Vec *v, MatNullSpace *nullSpace)
296: {
297: Vec vec;
298: PetscErrorCode (*funcs[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {zero_vector, constant_p};
301: DMGetGlobalVector(dm, &vec);
302: DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec);
303: VecNormalize(vec, NULL);
304: if (user->debug) {
305: PetscPrintf(PetscObjectComm((PetscObject)dm), "Pressure Null Space\n");
306: VecView(vec, PETSC_VIEWER_STDOUT_WORLD);
307: }
308: MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace);
309: if (v) {
310: DMCreateGlobalVector(dm, v);
311: VecCopy(vec, *v);
312: }
313: DMRestoreGlobalVector(dm, &vec);
314: /* New style for field null spaces */
315: {
316: PetscObject pressure;
317: MatNullSpace nullSpacePres;
319: DMGetField(dm, 1, &pressure);
320: MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullSpacePres);
321: PetscObjectCompose(pressure, "nullspace", (PetscObject)nullSpacePres);
322: MatNullSpaceDestroy(&nullSpacePres);
323: }
324: return 0;
325: }
327: int main(int argc, char **argv)
328: {
329: SNES snes; /* nonlinear solver */
330: DM dm, sdm; /* problem definition */
331: Vec u, r; /* solution, residual vectors */
332: Mat A, J; /* Jacobian matrix */
333: MatNullSpace nullSpace; /* May be necessary for pressure */
334: AppCtx user; /* user-defined work context */
335: PetscInt its; /* iterations for convergence */
336: PetscReal error = 0.0; /* L_2 error in the solution */
337: PetscReal ferrors[2];
338: PetscReal *coords, *viscosities;
339: PetscInt *materials;
340: const PetscInt particlesPerCell = 1;
341: PetscInt cStart, cEnd, c, d, bs;
344: PetscInitialize(&argc, &argv, NULL, help);
345: ProcessOptions(PETSC_COMM_WORLD, &user);
346: SNESCreate(PETSC_COMM_WORLD, &snes);
347: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
348: SNESSetDM(snes, dm);
349: DMSetApplicationContext(dm, &user);
351: PetscMalloc(2 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs);
352: SetupDiscretization(dm, &user);
353: //DMPlexCreateClosureIndex(dm, NULL);
355: /* Add a DMSWARM for particles */
356: DMCreate(PETSC_COMM_WORLD, &sdm);
357: DMSetType(sdm, DMSWARM);
358: DMSetDimension(sdm, user.dim);
359: DMSwarmSetCellDM(sdm, dm);
361: /* Setup particle information */
362: DMSwarmSetType(sdm, DMSWARM_PIC);
363: DMSwarmRegisterPetscDatatypeField(sdm, "material", 1, PETSC_INT);
364: DMSwarmRegisterPetscDatatypeField(sdm, "viscosity", 1, PETSC_REAL);
365: DMSwarmFinalizeFieldRegister(sdm);
367: /* Setup number of particles and coordinates */
368: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
369: DMSwarmSetLocalSizes(sdm, particlesPerCell * (cEnd - cStart), 4);
370: DMSwarmGetField(sdm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords);
371: DMSwarmGetField(sdm, "material", NULL, NULL, (void **)&materials);
372: DMSwarmGetField(sdm, "viscosity", NULL, NULL, (void **)&viscosities);
373: for (c = cStart; c < cEnd; ++c) {
374: const PetscInt i = (c - cStart) * bs;
375: PetscReal centroid[3];
377: DMPlexComputeCellGeometryFVM(dm, c, NULL, centroid, NULL);
378: for (d = 0; d < user.dim; ++d) coords[i + d] = centroid[d];
379: materials[c - cStart] = c % 4;
380: viscosities[c - cStart] = 1.0e20 + 1e18 * (cos(2 * PETSC_PI * centroid[0]) + 1.0);
381: }
382: DMSwarmRestoreField(sdm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords);
383: DMSwarmRestoreField(sdm, "material", NULL, NULL, (void **)&materials);
384: DMSwarmRestoreField(sdm, "viscosity", NULL, NULL, (void **)&viscosities);
386: DMCreateGlobalVector(dm, &u);
387: VecDuplicate(u, &r);
389: DMSetMatType(dm, MATAIJ);
390: DMCreateMatrix(dm, &J);
391: A = J;
392: CreatePressureNullSpace(dm, &user, NULL, &nullSpace);
393: MatSetNullSpace(A, nullSpace);
395: DMPlexSetSNESLocalFEM(dm, &user, &user, &user);
396: SNESSetJacobian(snes, A, J, NULL, NULL);
398: SNESSetFromOptions(snes);
400: DMProjectFunction(dm, 0.0, user.exactFuncs, NULL, INSERT_ALL_VALUES, u);
401: if (user.showInitial) DMVecViewLocal(dm, u, PETSC_VIEWER_STDOUT_SELF);
402: if (user.runType == RUN_FULL) {
403: PetscErrorCode (*initialGuess[2])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx) = {zero_vector, zero_scalar};
405: DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u);
406: if (user.debug) {
407: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
408: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
409: }
410: SNESSolve(snes, NULL, u);
411: SNESGetIterationNumber(snes, &its);
412: PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its);
413: DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);
414: DMComputeL2FieldDiff(dm, 0.0, user.exactFuncs, NULL, u, ferrors);
415: PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %.3g [%.3g, %.3g]\n", (double)error, (double)ferrors[0], (double)ferrors[1]);
416: if (user.showError) {
417: Vec r;
418: DMGetGlobalVector(dm, &r);
419: DMProjectFunction(dm, 0.0, user.exactFuncs, NULL, INSERT_ALL_VALUES, r);
420: VecAXPY(r, -1.0, u);
421: PetscPrintf(PETSC_COMM_WORLD, "Solution Error\n");
422: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
423: DMRestoreGlobalVector(dm, &r);
424: }
425: if (user.showSolution) {
426: PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
427: VecChop(u, 3.0e-9);
428: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
429: }
430: } else {
431: PetscReal res = 0.0;
433: /* Check discretization error */
434: PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
435: VecView(u, PETSC_VIEWER_STDOUT_WORLD);
436: DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error);
437: if (error >= 1.0e-11) PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", (double)error);
438: else PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n");
439: /* Check residual */
440: SNESComputeFunction(snes, u, r);
441: PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n");
442: VecChop(r, 1.0e-10);
443: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
444: VecNorm(r, NORM_2, &res);
445: PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", (double)res);
446: /* Check Jacobian */
447: {
448: Vec b;
449: PetscBool isNull;
451: SNESComputeJacobian(snes, u, A, A);
452: MatNullSpaceTest(nullSpace, J, &isNull);
454: VecDuplicate(u, &b);
455: VecSet(r, 0.0);
456: SNESComputeFunction(snes, r, b);
457: MatMult(A, u, r);
458: VecAXPY(r, 1.0, b);
459: VecDestroy(&b);
460: PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n");
461: VecChop(r, 1.0e-10);
462: VecView(r, PETSC_VIEWER_STDOUT_WORLD);
463: VecNorm(r, NORM_2, &res);
464: PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", (double)res);
465: }
466: }
467: VecViewFromOptions(u, NULL, "-sol_vec_view");
469: /* Move particles */
470: {
471: DM vdm;
472: IS vis;
473: Vec vel, locvel, pvel;
474: PetscReal dt = 0.01;
475: PetscInt vf[1] = {0};
476: PetscInt dim = user.dim, numSteps = 30, tn;
478: DMViewFromOptions(sdm, NULL, "-part_dm_view");
479: DMCreateSubDM(dm, 1, vf, &vis, &vdm);
480: VecGetSubVector(u, vis, &vel);
481: DMGetLocalVector(dm, &locvel);
482: DMGlobalToLocalBegin(dm, vel, INSERT_VALUES, locvel);
483: DMGlobalToLocalEnd(dm, vel, INSERT_VALUES, locvel);
484: VecRestoreSubVector(u, vis, &vel);
485: for (tn = 0; tn < numSteps; ++tn) {
486: DMInterpolationInfo ictx;
487: const PetscScalar *pv;
488: PetscReal *coords;
489: PetscInt Np, p, d;
491: DMSwarmVectorDefineField(sdm, DMSwarmPICField_coor);
492: DMCreateGlobalVector(sdm, &pvel);
493: DMSwarmGetLocalSize(sdm, &Np);
494: PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Timestep: %" PetscInt_FMT " Np: %" PetscInt_FMT "\n", tn, Np);
495: PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL);
496: /* Interpolate velocity */
497: DMInterpolationCreate(PETSC_COMM_SELF, &ictx);
498: DMInterpolationSetDim(ictx, dim);
499: DMInterpolationSetDof(ictx, dim);
500: DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
501: DMInterpolationAddPoints(ictx, Np, coords);
502: DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
503: DMInterpolationSetUp(ictx, vdm, PETSC_FALSE, PETSC_FALSE);
504: DMInterpolationEvaluate(ictx, vdm, locvel, pvel);
505: DMInterpolationDestroy(&ictx);
506: /* Push particles */
507: DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
508: VecViewFromOptions(pvel, NULL, "-vel_view");
509: VecGetArrayRead(pvel, &pv);
510: for (p = 0; p < Np; ++p) {
511: for (d = 0; d < dim; ++d) coords[p * dim + d] += dt * pv[p * dim + d];
512: }
513: DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
514: /* Migrate particles */
515: DMSwarmMigrate(sdm, PETSC_TRUE);
516: DMViewFromOptions(sdm, NULL, "-part_dm_view");
517: VecDestroy(&pvel);
518: }
519: DMRestoreLocalVector(dm, &locvel);
520: ISDestroy(&vis);
521: DMDestroy(&vdm);
522: }
524: MatNullSpaceDestroy(&nullSpace);
525: if (A != J) MatDestroy(&A);
526: MatDestroy(&J);
527: VecDestroy(&u);
528: VecDestroy(&r);
529: SNESDestroy(&snes);
530: DMDestroy(&sdm);
531: DMDestroy(&dm);
532: PetscFree(user.exactFuncs);
533: PetscFinalize();
534: return 0;
535: }