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, &user, &user, &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: }