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;
182: PetscFunctionBeginUser;
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: PetscCall(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: PetscCall(PetscOptionsEList("-bc_type", "Type of boundary condition", "ex62.c", bcTypes, 2, bcTypes[options->bcType], &bc, NULL));
195: options->bcType = (BCType)bc;
197: PetscCall(PetscOptionsBool("-show_initial", "Output the initial guess for verification", "ex62.c", options->showInitial, &options->showInitial, NULL));
198: PetscCall(PetscOptionsBool("-show_solution", "Output the solution for verification", "ex62.c", options->showSolution, &options->showSolution, NULL));
199: PetscCall(PetscOptionsBool("-show_error", "Output the error for verification", "ex62.c", options->showError, &options->showError, NULL));
200: PetscOptionsEnd();
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: PetscErrorCode DMVecViewLocal(DM dm, Vec v, PetscViewer viewer)
205: {
206: Vec lv;
207: PetscInt p;
208: PetscMPIInt rank, size;
210: PetscFunctionBeginUser;
211: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
212: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
213: PetscCall(DMGetLocalVector(dm, &lv));
214: PetscCall(DMGlobalToLocalBegin(dm, v, INSERT_VALUES, lv));
215: PetscCall(DMGlobalToLocalEnd(dm, v, INSERT_VALUES, lv));
216: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Local function\n"));
217: for (p = 0; p < size; ++p) {
218: if (p == rank) PetscCall(VecView(lv, PETSC_VIEWER_STDOUT_SELF));
219: PetscCall(PetscBarrier((PetscObject)dm));
220: }
221: PetscCall(DMRestoreLocalVector(dm, &lv));
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
226: {
227: PetscFunctionBeginUser;
228: PetscCall(DMCreate(comm, dm));
229: PetscCall(DMSetType(*dm, DMPLEX));
230: PetscCall(DMSetFromOptions(*dm));
231: PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view"));
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: PetscErrorCode SetupProblem(DM dm, AppCtx *user)
236: {
237: PetscDS prob;
239: PetscFunctionBeginUser;
240: PetscCall(DMGetDS(dm, &prob));
241: PetscCall(PetscDSSetResidual(prob, 0, f0_u, f1_u));
242: PetscCall(PetscDSSetResidual(prob, 1, f0_p, f1_p));
243: PetscCall(PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu));
244: PetscCall(PetscDSSetJacobian(prob, 0, 1, NULL, NULL, g2_up, NULL));
245: PetscCall(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: PetscFunctionReturn(PETSC_SUCCESS);
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;
271: PetscFunctionBeginUser;
272: /* Create finite element */
273: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
274: PetscCall(PetscFECreateDefault(comm, dim, dim, user->simplex, "vel_", PETSC_DEFAULT, &fe[0]));
275: PetscCall(PetscObjectSetName((PetscObject)fe[0], "velocity"));
276: PetscCall(PetscFECreateDefault(comm, dim, 1, user->simplex, "pres_", PETSC_DEFAULT, &fe[1]));
277: PetscCall(PetscFECopyQuadrature(fe[0], fe[1]));
278: PetscCall(PetscObjectSetName((PetscObject)fe[1], "pressure"));
279: /* Set discretization and boundary conditions for each mesh */
280: while (cdm) {
281: PetscCall(DMGetDS(cdm, &ds));
282: PetscCall(PetscDSSetDiscretization(ds, 0, (PetscObject)fe[0]));
283: PetscCall(PetscDSSetDiscretization(ds, 1, (PetscObject)fe[1]));
284: PetscCall(SetupProblem(cdm, user));
285: if (user->bcType == NEUMANN) PetscCall(DMGetLabel(cdm, "boundary", &label));
286: else PetscCall(DMGetLabel(cdm, "marker", &label));
287: PetscCall(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: PetscCall(DMGetCoarseDM(cdm, &cdm));
289: }
290: PetscCall(PetscFEDestroy(&fe[0]));
291: PetscCall(PetscFEDestroy(&fe[1]));
292: PetscFunctionReturn(PETSC_SUCCESS);
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};
300: PetscFunctionBeginUser;
301: PetscCall(DMGetGlobalVector(dm, &vec));
302: PetscCall(DMProjectFunction(dm, 0.0, funcs, NULL, INSERT_ALL_VALUES, vec));
303: PetscCall(VecNormalize(vec, NULL));
304: if (user->debug) {
305: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)dm), "Pressure Null Space\n"));
306: PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_WORLD));
307: }
308: PetscCall(MatNullSpaceCreate(PetscObjectComm((PetscObject)dm), PETSC_FALSE, 1, &vec, nullSpace));
309: if (v) {
310: PetscCall(DMCreateGlobalVector(dm, v));
311: PetscCall(VecCopy(vec, *v));
312: }
313: PetscCall(DMRestoreGlobalVector(dm, &vec));
314: /* New style for field null spaces */
315: {
316: PetscObject pressure;
317: MatNullSpace nullSpacePres;
319: PetscCall(DMGetField(dm, 1, &pressure));
320: PetscCall(MatNullSpaceCreate(PetscObjectComm(pressure), PETSC_TRUE, 0, NULL, &nullSpacePres));
321: PetscCall(PetscObjectCompose(pressure, "nullspace", (PetscObject)nullSpacePres));
322: PetscCall(MatNullSpaceDestroy(&nullSpacePres));
323: }
324: PetscFunctionReturn(PETSC_SUCCESS);
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;
343: PetscFunctionBeginUser;
344: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
345: PetscCall(ProcessOptions(PETSC_COMM_WORLD, &user));
346: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
347: PetscCall(CreateMesh(PETSC_COMM_WORLD, &user, &dm));
348: PetscCall(SNESSetDM(snes, dm));
349: PetscCall(DMSetApplicationContext(dm, &user));
351: PetscCall(PetscMalloc(2 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.exactFuncs));
352: PetscCall(SetupDiscretization(dm, &user));
353: //PetscCall(DMPlexCreateClosureIndex(dm, NULL));
355: /* Add a DMSWARM for particles */
356: PetscCall(DMCreate(PETSC_COMM_WORLD, &sdm));
357: PetscCall(DMSetType(sdm, DMSWARM));
358: PetscCall(DMSetDimension(sdm, user.dim));
359: PetscCall(DMSwarmSetCellDM(sdm, dm));
361: /* Setup particle information */
362: PetscCall(DMSwarmSetType(sdm, DMSWARM_PIC));
363: PetscCall(DMSwarmRegisterPetscDatatypeField(sdm, "material", 1, PETSC_INT));
364: PetscCall(DMSwarmRegisterPetscDatatypeField(sdm, "viscosity", 1, PETSC_REAL));
365: PetscCall(DMSwarmFinalizeFieldRegister(sdm));
367: /* Setup number of particles and coordinates */
368: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
369: PetscCall(DMSwarmSetLocalSizes(sdm, particlesPerCell * (cEnd - cStart), 4));
370: PetscCall(DMSwarmGetField(sdm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords));
371: PetscCall(DMSwarmGetField(sdm, "material", NULL, NULL, (void **)&materials));
372: PetscCall(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: PetscCall(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: PetscCall(DMSwarmRestoreField(sdm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords));
383: PetscCall(DMSwarmRestoreField(sdm, "material", NULL, NULL, (void **)&materials));
384: PetscCall(DMSwarmRestoreField(sdm, "viscosity", NULL, NULL, (void **)&viscosities));
386: PetscCall(DMCreateGlobalVector(dm, &u));
387: PetscCall(VecDuplicate(u, &r));
389: PetscCall(DMSetMatType(dm, MATAIJ));
390: PetscCall(DMCreateMatrix(dm, &J));
391: A = J;
392: PetscCall(CreatePressureNullSpace(dm, &user, NULL, &nullSpace));
393: PetscCall(MatSetNullSpace(A, nullSpace));
395: PetscCall(DMPlexSetSNESLocalFEM(dm, PETSC_FALSE, &user));
396: PetscCall(SNESSetJacobian(snes, A, J, NULL, NULL));
398: PetscCall(SNESSetFromOptions(snes));
400: PetscCall(DMProjectFunction(dm, 0.0, user.exactFuncs, NULL, INSERT_ALL_VALUES, u));
401: if (user.showInitial) PetscCall(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: PetscCall(DMProjectFunction(dm, 0.0, initialGuess, NULL, INSERT_VALUES, u));
406: if (user.debug) {
407: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n"));
408: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));
409: }
410: PetscCall(SNESSolve(snes, NULL, u));
411: PetscCall(SNESGetIterationNumber(snes, &its));
412: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
413: PetscCall(DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error));
414: PetscCall(DMComputeL2FieldDiff(dm, 0.0, user.exactFuncs, NULL, u, ferrors));
415: PetscCall(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: PetscCall(DMGetGlobalVector(dm, &r));
419: PetscCall(DMProjectFunction(dm, 0.0, user.exactFuncs, NULL, INSERT_ALL_VALUES, r));
420: PetscCall(VecAXPY(r, -1.0, u));
421: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solution Error\n"));
422: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
423: PetscCall(DMRestoreGlobalVector(dm, &r));
424: }
425: if (user.showSolution) {
426: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solution\n"));
427: PetscCall(VecFilter(u, 3.0e-9));
428: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));
429: }
430: } else {
431: PetscReal res = 0.0;
433: /* Check discretization error */
434: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n"));
435: PetscCall(VecView(u, PETSC_VIEWER_STDOUT_WORLD));
436: PetscCall(DMComputeL2Diff(dm, 0.0, user.exactFuncs, NULL, u, &error));
437: if (error >= 1.0e-11) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", (double)error));
438: else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: < 1.0e-11\n"));
439: /* Check residual */
440: PetscCall(SNESComputeFunction(snes, u, r));
441: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial Residual\n"));
442: PetscCall(VecFilter(r, 1.0e-10));
443: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
444: PetscCall(VecNorm(r, NORM_2, &res));
445: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", (double)res));
446: /* Check Jacobian */
447: {
448: Vec b;
449: PetscBool isNull;
451: PetscCall(SNESComputeJacobian(snes, u, A, A));
452: PetscCall(MatNullSpaceTest(nullSpace, J, &isNull));
453: //PetscCheck(isNull,PETSC_COMM_WORLD, PETSC_ERR_PLIB, "The null space calculated for the system operator is invalid.");
454: PetscCall(VecDuplicate(u, &b));
455: PetscCall(VecSet(r, 0.0));
456: PetscCall(SNESComputeFunction(snes, r, b));
457: PetscCall(MatMult(A, u, r));
458: PetscCall(VecAXPY(r, 1.0, b));
459: PetscCall(VecDestroy(&b));
460: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Au - b = Au + F(0)\n"));
461: PetscCall(VecFilter(r, 1.0e-10));
462: PetscCall(VecView(r, PETSC_VIEWER_STDOUT_WORLD));
463: PetscCall(VecNorm(r, NORM_2, &res));
464: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Linear L_2 Residual: %g\n", (double)res));
465: }
466: }
467: PetscCall(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: PetscCall(DMViewFromOptions(sdm, NULL, "-part_dm_view"));
479: PetscCall(DMCreateSubDM(dm, 1, vf, &vis, &vdm));
480: PetscCall(VecGetSubVector(u, vis, &vel));
481: PetscCall(DMGetLocalVector(dm, &locvel));
482: PetscCall(DMGlobalToLocalBegin(dm, vel, INSERT_VALUES, locvel));
483: PetscCall(DMGlobalToLocalEnd(dm, vel, INSERT_VALUES, locvel));
484: PetscCall(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: PetscCall(DMSwarmVectorDefineField(sdm, DMSwarmPICField_coor));
492: PetscCall(DMCreateGlobalVector(sdm, &pvel));
493: PetscCall(DMSwarmGetLocalSize(sdm, &Np));
494: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "Timestep: %" PetscInt_FMT " Np: %" PetscInt_FMT "\n", tn, Np));
495: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, NULL));
496: /* Interpolate velocity */
497: PetscCall(DMInterpolationCreate(PETSC_COMM_SELF, &ictx));
498: PetscCall(DMInterpolationSetDim(ictx, dim));
499: PetscCall(DMInterpolationSetDof(ictx, dim));
500: PetscCall(DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
501: PetscCall(DMInterpolationAddPoints(ictx, Np, coords));
502: PetscCall(DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
503: PetscCall(DMInterpolationSetUp(ictx, vdm, PETSC_FALSE, PETSC_FALSE));
504: PetscCall(DMInterpolationEvaluate(ictx, vdm, locvel, pvel));
505: PetscCall(DMInterpolationDestroy(&ictx));
506: /* Push particles */
507: PetscCall(DMSwarmGetField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
508: PetscCall(VecViewFromOptions(pvel, NULL, "-vel_view"));
509: PetscCall(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: PetscCall(DMSwarmRestoreField(sdm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords));
514: /* Migrate particles */
515: PetscCall(DMSwarmMigrate(sdm, PETSC_TRUE));
516: PetscCall(DMViewFromOptions(sdm, NULL, "-part_dm_view"));
517: PetscCall(VecDestroy(&pvel));
518: }
519: PetscCall(DMRestoreLocalVector(dm, &locvel));
520: PetscCall(ISDestroy(&vis));
521: PetscCall(DMDestroy(&vdm));
522: }
524: PetscCall(MatNullSpaceDestroy(&nullSpace));
525: if (A != J) PetscCall(MatDestroy(&A));
526: PetscCall(MatDestroy(&J));
527: PetscCall(VecDestroy(&u));
528: PetscCall(VecDestroy(&r));
529: PetscCall(SNESDestroy(&snes));
530: PetscCall(DMDestroy(&sdm));
531: PetscCall(DMDestroy(&dm));
532: PetscCall(PetscFree(user.exactFuncs));
533: PetscCall(PetscFinalize());
534: return 0;
535: }